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ABSTRACT 

The two fundamental assumptions of the standard cosmological model — that the initial fluctuations are statistically isotropic and Gaussian — are 
rigorously tested using maps of the CMB anisotropy from the Planck satellite. The detailed results are based on studies of four independent esti- 
mates of the CMB that are compared to simulations using a fiducial ACDM model and incorporating essential aspects of the Planck measurement 
process. Deviations from isotropy have been found and demonstrated to be robust against component separation algorithm, mask and frequency 
dependence. Many of these anomalies were previously observed in the WMAP data, and are now confirmed at similar levels of significance (around 
3cr). However, we find little evidence for non-Gaussianity with the exception of a few statistical signatures that seem to be associated with specific 
anomalies. In particular, we find that the quadrupole-octopole alignment is also connected to a low observed variance of the CMB signal. The 
dipolar power asymmetry is now found to persist to much smaller angular scales, and can be described in the low-^ regime by a phenomenological 
dipole modulation model. Finally, it is plausible that some of these features may be reflected in the angular power spectrum of the data which 
shows a deficit of power on the same scales. Indeed, when the power spectra of two hemispheres defined by a preferred direction are considered 
separately, one shows evidence for a deficit in power, whilst its opposite contains oscillations between odd and even modes that may be related to 
the parity violation and phase correlations also detected in the data. Whilst these analyses represent a step forward in building an understanding of 
the anomalies, a satisfactory explanation based on physically motivated models is still lacking. 

Key words. Cosmology: observations, cosmic microwave background, isotropy, Gaussianity 



1. Introduction 



This paper, one of a set associ ated with the 2013 release of 
data from the Pterc/Q mission ( [Planck Collaboration I 2013 ) 
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1 Planck |http : //www . esa . int/Plan"ck| > is a project of the 
European Space Agency (ESA) with instruments provided by two sci- 
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describes a set of studies undertaken to determine the statistical 
properties of the cosmic microwave background (CMB). 

The standard cosmological model is well described by the 
Friedmann-Lemaitre-Robertson-Walker (FLRW) solution of the 
Einstein field equations. This model is characterized by a ho- 
mogeneous and isotropic metric and an expanding scale factor 
of the Universe. At very early times it is hypothesized that the 
universe went through a period of accelerated expansion, the 
so-called cosmological inflation, driven by a hypothetical scalar 
field, the inflaton. During inflation the universe behaves as a 
de Sitter space, providing the conditions in which some of the 
present properties of the universe can be realized and specifi- 
cally relaxing the problem of initial conditions. In particular, the 
seeds that gave rise to the present large-scale matter distribution 
via gravitational instability originated as quantum fluctuations 
of the inflaton about its vacuum state. These fluctuations in the 
inflaton produce energy perturbations which are distributed as a 
homogeneous and isotropic Gaussian random field. Linear the- 
ory relates those energy fluctuations to the CMB anisotropies, 
implying a distribution for the anisotropies very close to that of 
an isotropic Gaussian random field. 

The scope of this work is to use Planck data to test the 
Gaussianity and near isotropy of the CMB in intensity as ex- 
pected in the standard cosmology paradigm. Testing these fun- 
damental properties is crucial for the validation of the standard 
cosmological scenario, and has profound implications for our 
understanding of the physical nature of the Universe and the ini- 
tial conditions of structure formation. Moreover, the confirma- 
tion of the isotropic and Gaussian nature of the CMB is essential 
to justify the corresponding assumptions usually made in the es- 
timation of the CMB power spectra, and other quantities to be 
obtained from the Planck data. Conversely, the detection of sig- 
nificant deviations from these assumptions that are not consistent 
with known systematic effects or foreground residuals would ne- 
cessitate major revision of current methodological approaches 
for the derivation of the mission's many science results. 

Significant deviations from Gaussianity are expected from 
non-linear processes that lead to secondary anisotropies, e.g. the 
integrated Sachs- Wolfe (ISW) effect and lensing. Indee d, these 
effects are the subject of two companion Planck papers ([Pl anck 



Collaboration XVII pOBl jPlanck Collaboratio n XIX 2qT3f7e 
spectively). However, remarkably, a number of anomalies, by 
which we mean features of the observed sky that are not statisti- 
cally consistent with the best-fit ACDM model, have been found 
in the WMAP data. Indeed, the WMAP team ( [Spergel et al.|2QQ3| ) 
themselves initially proposed some intriguing discrepancies in 
the form of a lack of power on large angular scales. Further ex- 
amples include an alignment of the low order multipoles some of 
which also indicate anomalously low amplitudes (Tegmark et al 



20031 |Bielewicz et al.|2005| |Land & Magueij o|2005a|), a NortfT 
South asymm etry in both power spectra ( [Eriksen et a l. 2004a; 



Han sen et al.||2009|) and various measures of non-Gaussianity 
( [Eriksen et al.|2004c|[2005l|Rath et al.|2007a| ), parity asymme- 



try in the power spectrum corresponding to large angular scales 
( Kim & Naselsky 2010a ) and a region o f significant temperature 
decrement — the so-called Cold Spot ( [Vielva et al.|2004[ [Cruz 
|et al.|2005| ). 

Whilst WMAP have presented refutations of these anoma- 
lies, either by criticism of the robustness of the statistical meth- 



ods employed ( Bennett et al.|2011 ) or by associating them with 
systematic artefacts of the data processing that have been cor- 
rected in the nine-year data release (Benne ttet al.|2012) , Planck 
represents a unique opportunity to independently assess their ex- 
istence. Its higher angular resolution and sensitivity and wider 
frequency range will allow a better understanding and removal 
of the Galactic and extragalactic foregrounds thus allowing a 
larger fraction of the sky to be useful for performing isotropy and 
Gaussianity analysis and to confirm and interpret those anoma- 
lies. 

Throughout this paper, we quantify the significance of the 
test statistic in terms of the p-value. This is the probability of 
obtaining a test statistic at least as extreme as the observed one, 
under the assumption that the null hypothesis (i.e., Gaussianity 
and isotropy of the CMB) is true. In some tests, where it is well 
motivated to use only a one-tailed probability, the p-value is re- 
placed by the corresponding upper or lower-tail probability. A 
low p-value is indicative of a tension between the data and the 
assumed statistical model (the null hypothesis). This can arise 
either when the assumed cosmological model is incorrect, if un- 
known or unmodelled aspects of the foreground emission or the 
measurement process exist, or as a result of a natural statistical 
fluctuation. The most interesting possibility, of course, is that a 
low p-value is an indication of new physics. 

From the theoretical point of view, there are many vari- 
ants of inflation that predict high levels of non-Gaussianity and 
new scenarios motivated by string and M-theory. In addition, 
there are many physical effects that might give rise to a devi- 
ation from isotropy or the presence of non-Gaussianity. Those 
deviations may be classified according to their physical na- 
ture and origin as follows: non-standard inflationary models, 
geometry and topology of the Universe, and topological de- 
fects. The main results from these areas, as well as the detailed 
descriptions of methodologies and of specific theoretically- 
motivate d model constraints, are provided in the co mpanio n pa- 
pers |Planck Collaboration XXIV|p013|), jPlanck Collaboration! 



entific consortia funded by ESA member states (in particular the lead 
countries France and Italy), with contributions from NASA, (USA) and 
telescope reflectors provided by a collaboration between ESA and a sci- 
entific conosrtium led and funded by Denmark 



|XXVI| ( |2013| ), and |Planck Collaboration XXV| ( [201 3[ T 

This paper covers all relevant aspects related to the phe- 
nomenological study of the statistical isotropy and Gaussian na- 
ture of the CMB measured by the Planck satellite. It is organized 
as follows. Section[2]describes the Planck data used for the anal- 
yses. Section [3] explains the main characteristics of the simula- 
tions that constitute our reference set of Gaussian sky maps rep- 
resentative of the null hypothesis. In Sect. [4] the null hypothesis 
is tested with a number of standard tests that probe different as- 
pects of non-Gaussianity. The WMAP anomalies are revisited in 
the light of the Planck data in Sect. [5 In Sect.[6]the implications 
of the found deviations of the null hypothesis on Q and cos- 
mological parameters estimations are discussed. Finally, Sect. [7] 
provides the main conclusions of the paper. 

2. Data description 

In this paper, we utilise data from the Planck-2013 data release 
corresponding to the nominal period of the Planck mission. In 
part, this comprises sky maps at nine frequencies, with corre- 
sponding 'half-ring' maps that are generated by separating the 
data for a given pointing period into two halves, plus maps gen- 
erated from data within the first and second survey periods. This 
set of maps allow a variety of consistency checks to be made, 
together with estimates of the instrumental noise contributions 
to analyses and limits on time- varying systematic artefacts. Full 
details are provided in papers JPlanck Collaboration II| ( |2013| ); 



Planck Collaboration^ ([2013^ 
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Our main results are based on the CMB maps resulting from 
sophisticated component separation algorithms applied to the 
frequency maps, as detailed in Planck Collabor ation XII| ( [2Q13j ). 
The four methods — Commander-Ruler, NILC, SEVEM and 
SMICA — are used to generate estimates of the CMB sky with an 
effective angular resolution of around 7' or better, with accompa- 
nying symmetrised beam profiles, analysis masks, half-ring and 
survey maps. In general, the analyses presented here make use 
of a standardised common mask that merges those associated 
with the individual methods (this mask is listed in Table [T] as 
U73). This is a conservative approach and therefore, where ap- 
propriate, we manipulate the masks for use at lower resolution. 
Low resolution maps are required in some analyses and have 
been produced as follows. For resolutions Af s id e = 128-1024 the 
full resoluti on maps have bee n degraded using the ud_de grade 
HEALPix (Gorski 'et al.||20 05) routine. For degrading to even 
lower resolutions, A^de =16-64, a different procedure has been 
followed. Before degrading the maps to the final resolution using 
the ud_degrade routine as in the previous case, the full-resolution 
map is smoothed with a Gaussian kernel with a FWHM equal to 
three times the pixel size of the low resolution map that we want 
to produce. 

In Table Q] we list the different masks that have been used 
in the analyses described in this paper. These masks have been 
produced at full resolution (A^ s i de =2048) and are described in 
papers |Planck Co llaborati on XII (2013); Planck Collaboration 
XV ( |2013| ). The mask U73 is the most often used in this paper. 
However, for several applications the masks have been degraded 
to lower resolutions (Af side =1024, 512, 256, 128, 64, 32 and 16). 
The masks with resolutions N S ^ Q = 128-1024 have been degraded 
using the following procedure: first, the masks are degraded to 
their final resolution using the ud_degrade HEALPix routine, 
and then, a conservative approach is followed setting to zero any 
pixels with a value lower than 0.8. If the masks have to be de- 
graded to even lower resolutions, A^de =16-64, the procedure 
that has been used is different. First, the full-resolution mask is 
smoothed with a Gaussian kernel with a FWHM equal to three 
times the pixel size of the low resolution mask that we want to 
produce. Then the mask is degraded using ud_degrade to their 
final resolution. Finally, those pixels with a value lower or equal 
than 0.5 have been set to zero and the rest have been set to 1 . This 
criterion, less conservative than the one used for the higher reso- 
lution masks, is a compromise between minimizing the amount 
of sky that is being masked and the level of contamination left 
unmasked (we remark that in some cases the more conservative 
criterion of a 0.8 threshold has been also used for the lower res- 
olutions, as stated in the corresponding analyses). 

3. Simulations 

The derivation of results to be presented in this paper requires ex- 
tensive simulations, essential aspects of which include: 1) mod- 
eling the Planck instrumental effects that affect the quality of 
the data, including instrumental noise and identified systematic 
effects, 2) replicating the foreground removal approach and es- 
timating the extent of foreground residuals, and 3) modeling 
the intrinsic statistical properties, Gaussian or otherwise, of the 
CMB signals expected from specific models of the Universe. 

The full focal plane (FFP6) simulations described in Planck 
Collaboration ES ( 2013 ) provide a complete realisation of the 
Planck mission capturing all characteristics of the scanning strat- 
egy, telescope, detector responses, and data reduction pipeline 
over the nominal mission period of 15.5 months. The Planck Sky 
Model (PSM) is used as input, encompassing the best current es- 



Table 1. List of the masks that have been used for the analy- 
ses described in this paper. All of them have been generated at 
A^ S ide = 2048, and when needed, they have been degraded to a 
lower resolution as explained in the text. The CL masks have 
been constructed following the procedure described in Planck 
Collaborati on XV| ( |2013| ) but for different sky coverages. 



Mask name 


Sky coverage 




[% of unmasked pixels] 



CS-SMICA89 1 

U73 1 

CG90 1 

CG80 1 

CG70 1 

CG60 1 

CL65 2 

CL58 2 

CL48 2 

CL37 2 

CL25 2 



89.0 
73.0 
90.0 
80.0 
70.0 
60.0 
65.1 
57.8 
48.0 
37.3 
24.7 



1 Planck Collaboration XII ( 2013 ) 
^Planck Coi laboratio nXV| ( |2013| ) 



timate of the microwave sky at Planck wavelengths including 
Galactic and extragalactic astrophysical foreground emission. 
The outputs include a complete set of maps for all detectors with 
accompanying half-ring and survey splits generated for a refer- 
ence CMB sky. These have been used to test and validate various 
analysis tools, employed in turn to evaluate the CMB component 
separation algorithms as applied to the data set. This also allows 
an FFP6-based estimate of the foreground residuals remaining 
in the CMB sky after component separation to be evaluated, and 
their impact on various statistical estimators quantified. 

An accompanying set of Monte Carlo simulations provides 
us with the reference set of Gaussian sky maps used f or the 
null tests we employ. These simulations include FEBeCoP (Mitra 
|et al.|[20TT] ) beam convolution at each of the Planck frequen- 
cies, which are then propagated through the various component 
separation pipelines using the same weights as derived from the 
Planck nominal mission data analysis. A fiducial CMB power 
spectrum has been adopted based on an analysis of the Planck 
data at an advanced, but not final stage of processing. Only small 
changes relative to the final Planck power spectrum presented 
in|Planck Collaboration XVl ( |20T3] ); [Pknck Collaboration XVI| 
( |2013| ) are observed. 



4. Are the primordial fluctuations Gaussian? 

As has been previously established, it is of major interest to de- 
termine whether the statistical properties of the primordial CMB 
anisotropies correspond to an isotropic Gaussian random field. 
Recent attempts to test this hypothesis have mainly relied on the 
WMAP data that have less sensitivity and cover a narrower fre- 
quency interval. Planck represents a unique opportunity to probe 
fundamental statistical properties of the Universe with cosmic 
variance limited sensitivity up to I « 2000 and minimum fore- 
ground contamination. 

There is no unique signature of non-Gaussianity, however, 
different tests can allow us to probe different types of non- 
Gaussianity. As a consequence, it is important to subject the 
data to a variety of tests, and we do so in this section using 
a number of non-parametric tools. Specific signatures of non- 
Gaussianity are sought in three companion papers — |Planck 



3 



Planck Collaboration: Isotropy and statistics 



Collaboratio n XXIV|([2QT3]);|Planck Collaboration XXVl ( |20T3] ); 
Planck Collaboration XXVllpQ13| ). 

Any isotropic and continuous random field, T(x) on the 
sphere can be written in terms of the following spectral repre- 
sentation: 



CO i 



T(X) = YjYj a ^m(x), 



(1) 



£=0 m =-l 



where x is a unit direction vector, Y{ m (-) the spherical harmonics 
and 



dim = J dxT(x)Y* m (x), 
m = 0,±l,...,±t, £ = 1,2,... 



For a Gaussian field with uncorrected phases, each a^ m coeffi- 
cient will be independent with a zero mean Gaussian distribu- 
tion: 



(at m a*, m ,) = 



(3) 



where S is the Kronecker delta and Q is the angular power spec- 
trum. Note that for a Gaussian and isotropic random field, the 
angular power spectrum provides a complete characterization of 
its statistical distribution. 

In this paper, we examine the goodness-of-fit of the data to 
the Planck best-fit fiducial CMB model, which constitutes our 
null hypothesis. The methods adopted constitute a broad range of 
statistical tools that allow the study of complementary statistical 
properties of the null hypothesis in both the real and harmonic 
space data representations. Claims of either consistency with 
the fiducial Planck cosmological model or of evidence for non- 
Gaussianity must be demonstrably robust to data selection and 
specifics of the data analysis. Residuals from the diffuse Galactic 
foreground are likely to be non-Gaussian in nature, and point- 
sources can be a source of non-Gaussianity on small angular 
scales. In addition, the analysis of multifrequency data must be 
considered in order to confirm that any claimed non-Gaussianity 
has a thermal (cosmological) spectrum. Moreover, the combined 
ISW-lensing effect produces secondary anisotropies that signif- 
icantly deviate from Gau ssianity and whose effect has be en de- 
tected in the Planck data ( [Planck Collaboration XIX|2Q13| ). This 
non-Gaussian effect has to be considered when testing the null 
hypothesis. 

We address these issues by analysing the cosmologically in- 
teresting subset of Planck frequency channels. Specifically, we 
analyse the uncorrected sky maps at 70, 100, 143 and 217 GHz 
as a function of Galactic mask to assess the likely contamina- 
tion due to Galactic foregrounds. These tests have dire ct rele- 
vance for the Plan ck likelihood approach described in Planck 
Collaboration XV|p013|) and the param eter estimation results 



presented in |Planck Col laboration XVI ( |2013 ). We then con 



sider the foreground cleaned versions of these maps generated 
by the SEVEM algorithm (see |Planck Collaboration XIl||20T3] ). 
Such a comparison also allows a semi-independent cross-check 
of the cosmological signal seen by Planck LFI (70 GHz) and 
HFI (100, 143, 217 GHz). Although the cosmological content of 
the cleaned LFI and HFI data sets are independent, the clean- 
ing makes use of difference maps generated from the remaining 
Planck frequency bands. Nevertheless, since the calibration and 
beam responses of the data are well understood over the ful range 
of frequencies, there will be no leakage of cosmological signal 
between the instrument specific frequencies. 

We then continue with analyses of the CMB sky es- 
timates provided by four component separation approaches 



Table 2. Lower tail probablity for the variance, skewness and 
kurtosis estimators at A^e = 2048, using the U73 mask and 
four different component separation methods. 



Method 


Variance Skewness 


Kurtosis 


C-R 


0.021 


0.189 


0.416 


NILC 


0.020 


0.191 


0.392 


SEVEM 


0.014 


0.206 


0.419 


SMICA 


. 0.017 


0.189 


0.419 



(Commander-Ruler, NILC, SEVEM, and SMICA) described in 
(2) |Planck C ollaboratio n XII| ( |2013| ), together with the corre- 



sponding mask appropriate for these methods. The largest sky 
area possible should be used for definitive statements about 
Gaussianity since, in the absence of foreground residuals or sys- 
tematic artefacts, it represents a superior sample of the Universe. 
Conversely, overly conservative sky cuts suffer from a loss of in- 
formation. 



4.1. One dimensional moments 

In this section we perform some of the simplest Gaussianity 
tests, such as comparing the sample skewness and kurtosis of 
the data with simulations. The skewness, y, and kurtosis, k, of a 
random variable, X, are defined as follows: 



7(X): 



K(X) 



(X-(X)) 3 
(Var(X)) 3 / 2 

(X-(X)) 4 
(Var(X)) 2 " 



(4) 



(5) 



The skewness is a measure of the asymmetry of the probabil- 
ity distribution of a real- valued random variable. Qualitatively, a 
positive (negative) skew indicates that the tail on the right (left) 
side of the probability density function is longer than the left 
(right) side. A zero value indicates that the values are relatively 
evenly distributed on both sides of the mean, typically but not 
necessarily implying a symmetric distribution. The kurtosis is a 
measure of the peakedness of the distribution and the heaviness 
of its tail. A distribution with positive (negative) excess kurto- 
sis indicates that the distribution has a more acute (wider) peak 
around the mean and fatter (thinner) tails. Normal random vari- 
ables have zero skewness and kurtosis. 

The sample variance is also considered in this section as a 
further consistency test, although it is not a normality test statis- 
tic. 

We begin by analysing the full resolution combined maps, 
applying the U73 mask for the four different component separa- 
tion methods. The results for the variance, skewness and kurtosis 
estimators are shown in Table |2] All four methods show similar 
results. The data are consistent with simulations for the skewness 
and kurtosis estimators, whereas the variance is anomalously 
lo w. This inconsistency was alre ady reported for th e WMAP data 
in |Monteserin et al.| ( |2008| ) and |Cruz et al.| ( |201l) at resolution 
Af side = 256 for a mask allowing slightly less sky coverage. 

The mask dependence of our results is studied by recalculat- 
ing the estimators using the CL58 and CL37 masks which allow 
sky fractions of / S k y = 58% and / S k y = 37% respectively. The 
SMICA cleaned maps at full resolution are considered. The most 
significant lower tail probability is obtained for the CL58 mask 
as can be seen in Table [3] The lower tail probabilities show a 
small dependence on the mask used, which could indicate ei- 
ther the presence of Galactic foreground residuals with larger 
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Table 3. Lower tail probablity for the variance, skewness and 
kurtosis estimators at A^de = 2048, for the SMICA method, using 
different masks. 



Mask 



Variance Skewness Kurtosis 



U73,/ sky =73% 

CL58, / sky =58% 

CL37, / sky =37% 

Ecliptic North, / sky =36% 
Ecliptic South, / sky =37% 
Galactic North, / sky =37% 
Galactic South, / sky =36% 



0.017 


0.189 


0.419 


0.003 


0.170 


0.363 


0.030 


0.314 


0.266 


0.001 


0.553 


0.413 


0.483 


0.077 


0.556 


0.001 


0.788 


0.177 


0.592 


0.145 


0.428 



Table 4. Lower tail probablity for the variance, skewness and 
kurtosis estimators at N S ^ Q = 2048, for the SEVEM cleaned maps 
at different frequencies. 



Map 


Variance Skewness Kurtosis 


100 GHz 


0.023 0.195 0.488 


143 GHz 


0.014 0.221 0.460 


217 GHz 


0.025 0.196 0.481 



sky coverage, or the increase of the sampling variance, and con- 
sequently a less significant probability, when a smaller fraction 
of the sky is considered. 

In order to identify any foreground contamination, the fre- 
quency dependence of our estimators is analysed. We use the 
SEVEM cleaned maps and the U73 mask. Note that as the 70 GHz 
full resolution noise is high we do not consider 70 GHz in this 
comparison. As the 100 GHz noise is not negligible we estimate 
th e variance taking in to account the noise dispersion as described 
in |Cruz e t al. (2011). The results are similar to those found for 
the combined map, as can be seen in Table [4] There is a small 
frequency dependence since the 100 GHz and 217 GHz maps 
show slightly higher variance and kurtosis than the 143 GHz 
map. However the 143 GHz map has a dominant contribution to 
the combined map, hence the foreground residuals in the com- 
bined map are likely to be small. The lower tail probabilities 
for the variance at 100 GHz, 143 GHz, and 217 GHZ are respec- 
tively 0.021, 0.014, 0.025, whereas the skewness and kurtosis 
are compatible with simulations. 

We also reanalyse the SMICA data and simulations consid- 
ering independently the northern and southern ecliptic hemi- 
spheres outside the U73 mask. A clear asymmetry is found in 
the variance, with an anomalously low value found in the north- 
ern hemisphere, as seen in Table [3] 

The results for different resolutions using the U73 mask are 
shown in Table[5] Note that the Af sid e = 2048 and 512 U73 masks 
have / S k y = 73%, while the low resolution masks at A^e = 64, 
32, and 16 have / S k y = 78%. The variance is anomalously low 
at all the considered resolutions, whereas at low resolutions, the 
skewness is anomalously low and the kurtosis anomalously high. 
These results will be further analysed in Sect. 5.2 However, it 
is clear that, except on the largest angular scales, there is no evi- 
dence for non-Gaussian behaviour in the data using these simple 
statistical measures. 



Table 5. Lower tail probablity for the variance, skewness and 
kurtosis estimators at different resolutions, for the four compo- 
nent separation methods, using the U73 mask. 



Method 


Variance Skewness Kurtosis 


Af side = 2048 


C-R 


0.021 


0.189 


0.416 


NILC 


0.020 


0.191 


0.392 


SEVEM 


0.014 


0.206 


0.419 


SMICA 


. 0.017 


0.189 


0.419 


Mide =512 


C-R 


0.017 


0.207 


0.368 


NILC 


0.017 


0.198 


0.390 


SEVEM 


0.013 


0.218 


0.408 


SMICA 


. 0.014 


0.196 


0.390 


Nside = 64 


C-R 


0.011 


0.041 


0.935 


NILC 


0.011 


0.041 


0.935 


SEVEM 


0.008 


0.058 


0.900 


SMICA 


. 0.011 


0.041 


0.943 


Wside = 32 


C-R 


0.020 


0.015 


0.968 


NILC 


0.019 


0.016 


0.960 


SEVEM 


0.012 


0.026 


0.932 


SMICA 


. 0.019 


0.016 


0.967 


Nside = 16 


C-R 


0.023 


0.013 


0.974 


NILC 


0.022 


0.022 


0.972 


SEVEM 


0.019 


0.022 


0.964 


SMICA 


. 0.027 


0.021 


0.982 



4.2. N-pdf analysis 

Under the assumption of Gaussianity, the Af -probability density 
function (N-pdf) is given by a multivariate Gaussian function: 



f(T) = 



1 



(27r) AW2 detC l/2 



exp-^rC" 1 ^), 



(6) 



where T is a vector formed from the measured temperatures T(x) 
over all positions allowed by the applied mask, N P i X is the num- 
ber of pixels in the vector, C is the covariance of the Gaussian 
field (of size N p [ x x N p [ x ). 

Unfortunately, the calculation of rC _1 r T is computationally 
unfeasible for the full Planck resolution at HEALPix N^q = 
2048. At a lower resolution, the problem is tractable, and the 
noise level can also be considered negligible compared to the 
CMB signal. That implies that under the assumption of isotropy 
the covariance matrix C is fully defined by the Planck angular 
power spectrum (Q): 



' max 



21 + 1 

47T 



QbjP, 



(cos%), 



(7) 



where Q 7 - is the covariance between pixels i and j, and 0^- is 
angle between them, Pi are the Legendre polynomials, hi is an 
effective window function associated with the Nside resolution, 
and ^ max is the maximum multipole probed. 

Under the multivariate Gaussian hypothesis, the argument 
on the exponential in equation [6] should follow a^ 2 distribution 
with A/pi x degrees of freedom, or, equivalently (for A/pi x » 1) a 



VVltll i v pix V^gl^O wi 1 1 

normal distribution N (N V] 
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Fig. 1. Variance, skewness and kurtosis for the combined map of 
the four different component separation methods. From top row 
to bottom row C-R, NILC, SEVEM, SMICA. 

Table 6. Lower tail probablity for the Af-pdf, using different 
masks. 



U73 



9200 9400 9600 9800 10000 10200 




Ecliptic North 



6600 6800 7000 7200 



4500 



4800 



5100 



CL37 



Ecliptic South 



4200 



4400 

x 2 



4600 



4500 



4800 

x 2 



5100 



Fig. 2. Af-pdf x 1 for the U73 mask, CL58, CL37, ecliptic North 
and ecliptic South. The different colours represent the four com- 
ponent separation methods, namely C-R (green), NILC (blue), 
SEVEM (red), and SMICA (orange). 

Table 7. Frequency dependence of the lower tail probablity for 
the Af-pdf, using different masks. 



Mask C-R NILC SEVEM SMICA 

U73, / sky =78% 0.027 0.028 0.019 0.030 

CL58,/ sky =58% 0.137 0.137 0.147 0.146 

CL37,/ Sky =37% 0.409 0.415 0.420 0.436 

Ecliptic North, / sky =39% 0.024 0.022 0.021 0.021 

Ecliptic South, / sky =39% 0.170 0.196 0.183 0.193 



Mask 



70 GHz 100 GHz 143 GHz 217 GHz 



U73, / sky =78% . 
CL58, / sky =58% 
CL37, / sky =37% 



Ecliptic South, f s 



sky 





0.037 


0.058 


0.013 


0.124 




0.169 


0.123 


0.111 


0.169 




0.422 


0.366 


0.376 


0.386 


=39% . 


0.028 


0.050 


0.015 


0.083 


=39% . 


0.225 


0.233 


0.166 


0.330 



We begin by analysing the x 2 quantity for low resolution 
maps at A^de = 32 and filtering with a 5° FWHM Gaussian. 
1 jjK uncorrected regularization noise is added to the co variance 
matrix before inverting it. Regularization noise realizations are 



added to the data and simulations for consistency (see Eriksen 
|et al.|2007b] for more details). 

We analyse the four cleaned data maps, applying the com- 
mon, CL58 and CL37 masks. The admitted fraction of the sky 
is respectively 78%, 58% and 37%. The northern and southern 
ecliptic hemispheres outside the U73 mask are also considered. 
The results are shown in Fig. [2] and Table [6] In the U73 mask 
case, the lower tail probabilities are low. Applying the two CL58 
and CL37 masks that permit less sky coverage, the data are con- 
sistent with simulations. The low^ 2 value appears to be localised 
in the northern ecliptic hemisphere. These results are directly 



comparable to the anomalous variance mentioned in Sect. 4.1 



Note that the four maps show similar values, but the differences 
are larger when using the U73 mask. This could indicate the 
presence of some residual foreground contamination near the 
Galactic plane. Therefore, the frequency dependence of our es- 
timator is analysed in order to identify any possible foreground 
contamination. The results are shown in Fig. [3] and Table [7] A 
moderate frequency dependence is found when using the U73 
mask, which could indicate the presence of some foreground 



residuals near the Galactic plane. The frequency dependence of 
the results vanishes when using the CL58 and CL37 masks that 
exclude more of the sky from analysis. 

4.3. N -point correlation functions 

In this section we present tests of the non-Gaussianity of the 
Planck CMB maps using real-space Af-point correlation func- 
tions. While harmonic- space methods are often preferred over 
real- space methods for studying primordial fluctuations, real- 
space methods may have an advantage with respect to system- 
atics and foregrounds, since such effects are usually localized in 
real space. It is therefore important to analyse the data in both 
spaces in order to highlight different features. 

An Af-point correlation function is by definition the average 
product of N temperatures, measured in a fixed relative orienta- 
tion on the sky, 



C N (0 l9 . . . , 62N-3) = ATXnO • • • Ar(n^) 



(8) 



where the unit vectors fii, . . . , n# span an Af-point polygon on 
the sky. By assuming statistical isotropy, the Af-point functions 
are only functions of the shape and size of the Af-point poly- 
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Fig. 3. Frequency dependence for 70 GHz (green), 100 GHz 
(blue), 143 GHz (red) and 217 GHz (orange), and different 
masks. 

gon, and not on its particular position or orientation on the sky. 
Hence, the smallest number of parameters that uniquely deter- 
mines the shape and size of the Af -point polygon is 2N - 3. In 
practice, the functions are estimated by simple product averages 
over all sets of N pixels fulfilling the geometric requirements set 
by 0i , ... , #2v-3 characterising the shape and size of the polygon 



C#(0i,...,02Ar_ 3 ) = 



Y Ji (w[---w i N )(AT[---Ar N ) 



(9) 



Pixel weights w\, • • • , w l N can be introduced in order to reduce 
noise or mask boundary effects. Here they represent masking by 
being set to 1 for included pixels and to for excluded pixels. 

The main difficulty with computing Af-point functions is 
their computational scaling. The number of independent pixel 
combinations scales as (9(A^ x ), and for each combination of 
N pixels, 2N - 3 angular distances must be computed to 
uniquely determine the properties of the corresponding polygon. 
Computing the full Af -point function for N > 2 and N p { x > 10 5 is 
therefore computationally challenging. However, it is not neces- 
sary to include all possible Af-point configurations in order to 
produce interesting results. For instance, one may focus only 
on small angular scales, or on configurations with some spe- 
cial symmetry prope rties. By using the methods described by 
Erik sen et al.| ( |2004b] ), the computational expense then becomes 
tractable, since no CPU time is spent on excluded configurations. 
In this paper several such subsets are computed, covering three 
distinct ranges of scales, namely small (up to 3°), intermediate 
(up to 10°) and large angular scales (the full range between 0° 
and 180°). The shapes of considered polygons selected for the 
analysis are the pseudo-collapsed and equilateral configurations 
for the 3 -point function, and the the rhombic configuration com- 
posed of two equilateral triangles for the 4-point function. In the 
following, all results refer to the reduced 4-point function, i.e., 
corrected for the Gaussian contribution due to the Wick's theo- 
rem. The size of the polygons is parametrised by the length of 



the longer side of the triangle in the case of the pseudo-collapsed 
configuration, and the length of the side for the equilateral trian- 
gle and rhombus. 

We analyse the CMB estimates downgraded to N s ^q = 64 
and A^ S ide = 512 as well as at the original resolution of Af s ide = 
2048. In the case of the analysis at N^q = 64 the maps were 
additionally smoothed with FWHM of 165' (three times the pixel 
size for the downgraded map). Due to computational limitations, 
an analysis is possible on the full sky only in the case of reso- 
lution A^ S ide = 64. For the higher resolution maps, we perform 
the analysis on a set of non-overlapping discs. For A^ s id e = 512 
we uniformly retain, after masking, part of the sky with approx- 
imately 100 discs of radius 10°. Analogously to the analysis by 
Eriksen et al. (2005}, we consider two disc sets A and B with a 
relative offset between their grids such that the centres of the 
discs of set B are located in parts of the sky not covered by 
disc set A (see Fig. [4]). For studies at the original resolution 
A^side = 2048, we restrict the analysis to 20 discs with a radius of 
3° located randomly on an unmasked part of the sky (Fig. [4}. 




Fig. 4. Two sets of discs, A and B, each of radius 10° for the 
A^side = 512 CMB estimates (upper and middle figure, respec- 
tively) and a set of 20 randomly placed discs of radius 3° super- 
imposed on the U73 mask (blue region) for the CMB estimates 
at Af sid e = 2048 (lower figure). 



7 



Planck Collaboration: Isotropy and statistics 



As in Eriksen et al. ([2005 ), we consider the AZ-point correla- 
tion functions averaged over the disc sets. In order to minimize 
correlations between the discs, we subtract from the maps at res- 
olutions A^side = 512 and A^de = 2048 the best-fit multipoles 
computed for the ranges i < 18 and i < 60, respectively. This 
procedure corresponds in practice to a high-pass filtering of the 
maps. 

The low resolution versions of the U73 mask described ear- 
lier were used as required. Residual monopole and dipole con- 
tributions were then removed from the maps. 

A simple x 2 test is chosen to quantify the degree of agree- 
ment between the simulations and the observations, where x 2 as 
usual is defined by 



N hiD 



X 



£ (c N (6d - (c N m)) m;/ (c N (6j) - (c N (6j))) 
i,j=l 



(10) 



Here C^iOi) is the Af-point correlation function for i - th bin of 
separation angle, ^CaK#/)) is the corresponding average from 
the Monte Carlo (MC) ensemble, and 



{c N (6d)) (cffWj) - (c N (0j))) (11) 



is the covariance matrix. Although the inverse of the covariance 
matrix constructed from MC simulations can be biased, it is rel- 
atively small for 1000 simulations and has a negligible impact 
on the significance levels estimated from the simulations, as de- 
scribed below. 

This statistic is optimized for studying Gaussian distributed 
data. However, usually it also works quite well for mildly 
non-Gaussian distributions, and in particular symmetric ones. 
Nevertheless, as for any statistic constructed from MC simula- 
tions, it can also be used for non-Gaussian and asymmetrically 
distributed data. Below, we quote the significance level in terms 
of the fraction of simulations with a larger^ 2 value than the ob- 
served map. 

We analyse the mask dependence of the non-Gaussianity of 
the maps using the pseudo-collapsed 3-point correlation func- 
tion. The function averaged over disc set A is shown in Fig. [5] 
The corresponding probabilities of obtaining values of the 
statistic for the concordance ACDM model at least as large as 
the observed values are given in Table [5] 

Table 8. Probabilities of obtaining values of the x 2 statistic for 
the concordance ACDM model at least as large as the observed 
values of the statistic for the raw 143 GHz (first row) and fore- 
ground corrected 143 GHz SEVEM CMB maps (second row). 
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Fig. 5. The pseudo-collapsed 3-point function averaged over disc 
set A for the raw (upper figure) and SEVEM foreground corrected 
(lower figure) 143 GHz map at Af s ide = 512. Estimates of the 
multipoles for £ < 18 are removed from the sky maps. The black 
solid line indicates the mean for 1000 MC simulations and the 
shaded dark and light grey regions indicate the 68% and 95% 
confidence regions, respectively, for the CG90 (/ S k y = 0.9) mask. 

Table 9. Probabilities of obtaining values for the x 2 statistic of 
the A/-point functions shown in Fig. [6] for the Planck fiducial 
ACDM model at least as large as the observed values of the 
statistic for the Planck CMB maps with resolution parameter 
Nside = 64 estimated using the C-R, NILC, SEVEM and SMICA 
methods. 





C-R 


NILC 


SEVEM 


SMICA 


2-pt 


0.883 


0.859 


0.884 


0.855 


pseudo-coll. 3-pt. . . 


. . 0.922 


0.918 


0.945 


0.908 


equil. 3-pt 


0.962 


0.966 


0.978 


0.968 


4-pt 


. . 0.975 


0.977 


0.979 


0.977 



In summary, the pseudo-collapsed 3-point function does not 
show any significant deviation from Gaussianity for the raw 
143 GHz map masked with the CG60 (/ sky = 0.6) and CG70 
(/sky = 0.7) masks. To a lesser extent, this is true also for the 
CG80 (/ S k y = 0.8) mask. We do not see any significant deviation 
for any of the analysed masks after cleaning the 143 GHz map 
using the SEVEM method. 

The correlation functions for the four component separation 
methods and resolution parameters A^ s id e = 64, A^ s id e = 512 and 



A^ S i d e = 2048 are shown in Fig.[6} Fig.|7](disc set A), Fig.[8](disc 
set B) and Fig. [9j respectively. The probabilities of obtaining 
values of the^ 2 statistic for the Planck fiducial ACDM model at 
least as large as the observed values are given in the Tables [9] [T0| 
and[TT] respectively. 

The results show consistency between the CMB maps es- 
timated using the different component separation methods. We 
did not find statistically significant deviations of the CMB maps 



8 



Planck Collaboration: Isotropy and statistics 




20 40 60 80 100 120 20 40 60 

6 [deg] 9 [deg] 



Fig. 6. The 2-point (upper left), pseudo-collapsed (upper right), equilateral (lower left) 3-point and reduced rhombic 4-point (lower 
right) functions for the A^de = 64 CMB estimates. The black solid line indicates the mean for 1000 MC simulations and the shaded 
dark and light grey regions indicate the 68% and 95% confidence regions, respectively. 



Table 10. Probabilities of obtaining values of the x 2 statistic of 
the Af -point functions shown in Figs. [7] and [8] for the Planck fidu- 
cial ACDM model at least as large as the observed values of 
the statistic for Planck CMB maps with resolution parameter 
Nside = 512 estimated using the C-R, NILC, SEVEM and SMICA 
methods. 



C-R NILC 


SEVEM 


SMICA 


Two-point function 


A set 0.858 0.902 

Bset 0.351 0.370 


0.886 
0.404 


0.904 
0.376 


Pseudo-collapsed three-point function 


A set 0.568 0.565 

Bset 0.483 0.526 


0.651 
0.550 


0.603 
0.540 


Equilateral three-point function 


A set 0.004 0.032 

Bset 0.452 0.485 


0.045 
0.443 


0.043 
0.479 


Rhombic four-point function 


A set 0.104 0.102 

Bset 0.521 0.569 


0.102 
0.537 


0.107 
0.579 



Table 11. Probabilities of obtaining values of the x 2 statistic 
of the N-point functions shown in Fig. [9] for the Planck fidu- 
cial ACDM model at least as large as the observed values of 
the statistic for Planck CMB maps with resolution parameter 
Af side = 2048 estimated using the C-R, NILC, SEVEM, and SMICA 
methods. 





C-R 


NILC 


SEVEM 


SMICA 


2-pt 


0.335 


0.474 


0.573 


0.497 


pseudo-coll. 3-pt. . . 


. . 0.522 


0.463 


0.469 


0.448 


equil. 3-pt 


0.853 


0.789 


0.819 


0.796 


4-pt 


. . 0.532 


0.534 


0.579 


0.526 



from Gaussianity for any of the analysed scales. However, it is 
clear that the CMB maps smoothed and downgraded to N S ^ Q = 
64 show the largest deviation, especially for the 4-point correla- 
tion function, in comparison to the intermediate and small angu- 
lar scale analyses. 

4.4. Minkowski fu actionals 

Minkowski functionals (Minkowski 1903, hereafter MFs) de- 
scribe the morphology of fields in any dimension and have long 
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Fig. 7. The 2-point (upper left), pseudo-collapsed (upper right), equilateral (lower left) 3-point and reduced rhombic 4-point (lower 
right) functions averaged over disc set A for the A^de = 512 CMB estimates. Estimates of the multipoles for £ < 18 are removed 
from the sky maps. The black solid line indicates the mean for 1000 MC simulations and the shaded dark and light grey regions 
indicate the 68% and 95% confidence regions, respectively. 



been used as estimators of non-Gaussianity and anisotropy in the 
CMB (see e.g. |Gott et al.|1990||Mecke et al.|1994l|Schmarzing 
& Gorski|1998HKomatsu et al.|2003[|Eriksen et al.|2004cl|C urto 
et al.|2007[|De Troia et al 12007 ||Spergel et al.|2007[|Curto et al 
2008} |Hikage et al.|2008||Komatsu et al.|2009j ). They are addi- 
tive for disjoint regions of the sky and invariant under rotations 
and translations. Traditionally in the literature, the contours are 
defined by a threshold v, usually given in units of the sky stan- 
dard deviation (cr ). We compute MFs for the regions colder 
and hotter than a given threshold v. Thus, the three MFs, the 
area Vq(v) = M y )> the perimeter Vi(v) = C(v) and the genus 
V2M = G(v), are defined respectively as: 



Vb(v) = A(v) = 



Vi(v) = C(v) = 



Nr 



pix 



4A„ 



V 2 (v) = G(v) = ^-^(^hot - A^coid), 



(12) 



(13) 



(14) 



where Af y is the number of pixels where AT/o-q > v, N V { X is 
the total number of available pixels, A tot is the total area of the 



available sky, A^hot is the number of compact hot spots, A^ co id is 
the number of compact cold spots and S t is the contour length of 
each hot spot. We construct a fourth functional V^iy) = Cluster (v) 
w hich corresponds to N co \^ for negative v and A^hot for positive 
y ( [Ducout et al. 2012 ). Analytical expressions for a Gaussian 
random field can be derived in terms of v (see e.g. Vanmarcke 
1 1 98 3 [ |Matsubara|20 1 0| ) and give the following, 

Vkiy) = A k v k {y\ (15) 



with 

vjfc(v) = exp(-y 2 /2)^_!(y), k<2 

-v 2 

v 3 (v) = 



erfc (y/ V5) ' 



and 



H n {y) = e 



v 2 /2 



dy 



-v 2 /2 



(16) 
(17) 

(18) 



The amplitude A k depends only on the shape of the power spec- 
trum Q: 



1 



C02 



(2tf)(*+D/2 CJ 2 - k CJ k \ V2CT 



k<2 



(19) 
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Fig. 8. As Fig.^for disc set B. 



e [deg] 



(20) 



where ojk = 7i k ^ 2 /T(k/2 + 1), which gives ojq = I, QJ\ = 2, 
a>2 = 7i and o"o and cr\ are respectively the rms of the field and 
its first derivatives. These analytical expressions represent use- 
ful descriptions of the MFs which, for the case of a Gaussian 
random field, can be factorized as a function of the threshold 
and another of the shape and amplitude of the Q. We will use 
both the unnormalized (Vk) and unormalized (v&) MFs in the 
Gaussianity tests performed in this section. The unnormalized 
functionals are comp uted with a code th at was used for the analy- 
sis of Archeops data ( |Curto et al.|2 007 ) and has been thoroughly 
validated with Planck simulations, while for the normalized ones 
a code adapted to the high resolution Planck data and described 
in |Ducout et al.| ( |20T2] ) is used. 

By combining the MFs curves into a vector y of size 
n = ^thresholds x Junctionals » a null hypothesis test can be performed 
using &X 2 statistic given by: 

x 2 (y) = (y- (yG)) T c- l (y - < jG » (21) 

where y represents the MFs of the data, Jg those of the simu- 
lations and C is the covariance matrix. In order to assure con- 
vergence, in the case of the four normalized MFs C is estimated 
from 10 4 Gaussian simulations, drawn from the Planck fiducial 
power spectrum, having the same instrumental properties of ef- 
fective beam and noise as the data, the same applied mask and 



which have been processed in the same way to reach the cor- 
responding resolution. For the three unnormalized MFs, C was 
estimated from only 10 3 FFP6 simulations that proved to be suf- 
ficient for convergence. We compare the /fcpianck obtained from 
the data to the x 2 obtained from those simulations, and report 
the probability of having a value of x 2 larger than the mea- 
sured one, P (x 2 > ^pi anck )- We explore different resolutions rep- 
resented by the parameter A^e, different methods of component 
separation (Commander-Ruler, NILC, SEVEM, and SMICA) and 
different sky coverages. 

First, the three unnormalized MFs (Vk as a function of v, 
k = 0, 1, 2) are used to construct a test of the null hypothesis. 
The test assesses not only the primordial Gaussian hypothesis, 
but also whether the data is correctly represented by the sim- 
ulations in terms of power spectrum, systematics and the lens- 
ing effect. A set of 17 thresholds between -4 and +4 in steps 
of 0.5 are considered. The comparison between the MFs of the 
data provided by the four component separation methods and 
those corresponding to each of the four sets of 10 3 FFP6 sim- 
ulations representing each method, for the standard U73 mask, 
are shown in Fig. 10 From that figure, a deviation at a level of 
« 2cr can be seen for the contour and genus curves at a reso- 
lution A^ S ide = 512. The situation is very similar for the analy- 
ses performed at other resolutions, A^de = 1024, 256 and 128. 
Although the deviation is not particularly compelling because 
of the correlations among neighbouring thresholds, it is worth 
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Fig. 9. The 2-point (upper left), pseudo-collapsed (upper right), equilateral (lower left) 3-point and reduced rhombic 4-point (lower 
right) functions averaged over the disc set for the A^de = 2048 CMB estimates. Estimates of the multipoles for £ < 60 are removed 
from the sky maps. The black solid line indicates the mean for 1000 MC simulations and the shaded dark and light grey regions 
indicate the 68% and 95% confidence regions, respectively. 



mentioning that a possible explanation is the background of un- 
resolved sources that has been detected in Planck data with th e 
bispectrum estimators (see Planck Collaboration XXIV (2013)). 
In order to understand the effect of unresolved sources on the 
MFs, we added the point source residuals derived from the FFP6 
simulations as processed by the SEVEM algorithm to 100 real- 
isations which were then analysed. We conclude that the back- 
ground of unresolved sources may be responsible for at least part 
of the excess signal that is detected. The corresponding probabil- 
ities P (x 2 > ^p lanck ) derived from the MF values for each of the 
four component separation methods and resolutions are given in 
Table [12] The full resolution maps have been degraded to the 
lower resolution ones following the procedure described in sec- 
tion^ All the cases considered are compatible with the null hy- 
pothesis. 

In the second case, the four normalized MFs Vk = Vk/Ak 
(k = 0, 1, 2, 3) are used for the null hypothesis test. A set of 
26 thresholds equally spaced between -3.5 and +3.5 are con- 
sidered. The normalization factor is estimated directly from 
the maps, having computed previously the moments cr and <j\ . 
This normalization minimizes the dependence of the MFs on the 
power spectrum, thereby decreasing the cosmic variance and im- 
proving their sensitivity to deviations from Gaussianity. The res- 



Table 12. Non-directional Gaussianity tests using unnormalized 
MFs: P (x 2 > Afpi anc k) as a function of sky resolution for different 
component separation methods. 



side 


1024 


512 


256 


128 


C-R 


0.812 


0.299 


0.482 


0.357 


NILC 


0.993 


0.567 


0.354 


0.234 


SEVEM 


0.925 


0.911 


0.738 


0.094 


SMICA 


0.874 


0.675 


0.426 


0.213 



olutions considered in this case are N S id e = 2048, 1024, 512, 256 
and 128. For the highest resolution N s ^ e = 2048, the map 
is smoothed with a Gaussian smoothing kernel with a width 
#fwhm = 5', in order to decrease the noise level. We use the 
standard U73 mask, inpainting the smallest point sources. The 
maps at lower resolution are constructed by the standard simple 
degrading process applied to the original map at A^e = 2048, 
and the corresponding masks are degraded following a conserva- 
tive procedure such that any degraded pixel with a value < 0.8 is 
set to zero (as explained in section [2]). The results of the analysis 
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Fig. 10. Difference of the data MFs (unnormalized) with respect to the average of the curves obtained with realistic Planck sim- 
ulations for several cleaned maps. From left to right: Area, Contour, Genus. The error-bars represent the l<x (68%CL) dispersions 
around the mean obtained with simulations. 



performed on the SMICA map at different resolutions are pre- 
sented in Table [13] The results of the analysis performed on the 



Table 13. Non-directional Gaussianity tests using normalized 
Minkowski Functionals: dependence of p(x 2 > Api anc k) on Sky 
resolution. 



A^side 


2048 


1024 


512 


256 


128 


Normalized MFs 


0.358 


0.356 


0.245 


0.225 


0.223 



different component separation methods at the highest resolution 
(A^side = 2048) are presented in Table 14 The difference of the 
normalized MFs with respect to the expected values of the null 



hypothesis as a function of the threshold v are shown in Fig. 1 1 
A slight deviation in N c \ ustQY (v) is noticeable at thresholds v « 
however it is not very compelling since the values at neighbor- 
ing thresholds are very correlated and this correlation is taken 
into account in the^ 2 statistics. Finally, we analyse the depen- 



Table 14. Non directional Gaussianity tests using normalized 
Minkowski Functionals: Dependence on omponent separation 
methods. 



Method 


C-R 


NILC 


SEVEM 


SMICA 


P{)( > Apianck) 


0.288 


0.303 


0.415 


0.358 



dence of the normalized MFs on the sky coverage. We use the 
standard U73 mask and then decrease the sky coverage by us- 
ing CL65, CL48 and CL25 masks in combination with a special 
point source mask that is based on the U73 mask. The fraction 
of sky left unmasked in the combined masks is 62%, 46% and 
23%, respectively. The point source mask is inpainted previously 
to the analysis. The curves obtained for the different sky cover- 
ages are presented in Fig. 12 for the SMICA method. Results of 
the^ 2 analysis of the data as a function of sky coverage are com- 
piled in Table [15] All the cases considered are compatible with 
the null hypothesis. 

In summary, we find that the data are globally consistent with 
the primordial Gaussian hypothesis, and no strong deviation is 



found between the data and realistic simulations for both the un- 
normalized and normalized MFs. We would like to remark that 
a certain level of non-Gaussianity is expected from lensing and, 
in particular, from the ISW-lensing signal, thus it is important to 
compare the data to realistic lensed simulations. 



Table 15. Non directional Gaussianity tests using normalized 
Minkowski Functionals : Sky coverage. 



/sky 


0.73 


0.62 


0.46 


0.23 


P{)( > Apianck) 


0.358 


0.042 


0.670 


0.780 



4.5. Wavelet statistics 

A broad range of wavelets have been used in the analysis of 
CMB data, but in this paper we consider the Spheric al Mexican 
Hat wavelet (SMHW, [Martinez-Gonzalez et al.|2002) . 

The SMHW is an example of a continuous, non-orthogonal 
wavelet. Given a signal on the sky, T(p), where p represents a 
given position/pixel wich is a function of the co-latitude and 
longitude (also defined by the unit direction vector x), the 
SMHW coefficients at a given scale R, oj t (R, p), are obtained 
by convolution: 



(=0 m=-C 



SMHW 



(R)Y (m (p), 



(22) 



where w^ MHW (R) is the window function associated with the 
SMHW, 

^max is the maximum multipole allowed by the cor- 
responding HEALPix pixelization, Y{ m (p) is the spherical har- 
monic basis, and fy m are the spherical harmonic coefficients of 
the analysed map: 



t In 



-I 



dQY* m (p)T(p), 



(23) 



where dQ = d6 sin 6d(f) and * denotes complex conjugation. 

Several statistics can be computed from the wavelet coeffi- 
cients map, in particular, the first moments: the dispersion cr R , 
the skewness Sr, and the kurtosis K R (as a function of scale R). 
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v v 

Fig. 11. Difference of the normalized MFs obtained from the data with respect to the expected values of the null hypothesis, for the 
different component separation methods. From left to right and top to bottom: Area, Contour, Genus and Af c i uste r. The grey bands 
represent the 1 and 2<x dispersions around zero, based on realistic Planck simulations including lensing, for C-R method. 



It is interesting to notice that in the case of Gaussian tempera- 
ture fluctuations the linear transformatio n involv ed in the deter- 
mination of the wavelet coefficients (eqs. 22|23| ) guarantees that 
Gaussianity is preserved. 

The study of the moments of the distribution of the CMB 
temperature fluctuations, as a function of the scale, is a standard 
approach to test the null hypothesis. We have performed a full 
resolution multi-scale analysis of the four CMB clean maps and 
computed the quantities cr R , and from the SMHW wavelet 
coefficients at 18 scales, R = {2, 4, 7, 14, 25, 50, 75, 100, 150, 
200, 250, 300, 400, 500, 600, 750, 900, 1050}, in arc-minutes. 
These are compared to the standard Planck simulations. 



As explained in Viel va et"aL] ( |2004| ), when computing the 
SMHW coefficients of a masked data set, artefacts are intro- 
duced close to the mask that degrade the performance of any null 
hypothesis tests. We therefore define a set of exclusion masks 
such that, at each scale, an extra region of the sky is excluded 
when performing any statistical test. The exclusion mask for a 
given scale R is defined as follows: we build an auxiliary mask 
by removing from the U73 mask all the features associated with 
compact objects, and degrade this auxiliary mask to A^de = 1024 
(imposing a restrictive cut); a first temporary mask is obtained 
by extending the borders of this auxiliary mask by a distance 
of twice R; a second temporary mask is obtained, first, by con- 
volving the auxiliary mask with the SMHW at that particular 
scale R and, second, by imposing that any pixel of that sec- 
ond temporary mask with an absolute value lower than 0.1 is 
masked, whereas the remaining ones are set to 1 ; the two tempo- 
rary masks are multiplied to yield a single mask that is upgraded 



to A^ S ide = 2048; finally, the final exclusion mask is obtained by 
multiplying this mask by the parent U73 mask. 

The comparison of the four CMB maps with the correspond- 
ing simulations is summarized in Fig. [13] The three panels show 
(from left to right) the statistical significance of the standard 
deviation, the skewness and the kurtosis (as a function of the 
SMHW scale). The points represent the upper tail probability 
associated to a given statistic, i.e., the fraction of the simulations 
that present a value of a given statistic equal to or greater than 
the one obtained for the data. In fact, we define a modified upper 
tail probability: if an upper tail probability is larger than 0.5, then 
a new quantity is defined as 1 minus that upper tail probability. 
Hence, this modified definition of the upper tail probability is 
constrained between 10" 3 (the minimum value that can be im- 
posed with 1 000 simulations) and 0.5. Overall, the agreement 
between the four CMB maps is quite good, showing that all of 
them provide a consistent estimation of the true CMB. However, 
several aspects need to be discussed. Let us clarify that the dif- 
ferences among the CMB methods for small modified upper tail 
probabilities are expected to be larger than for large modified 
upper tail probabilities. This is because a small modified upper 
tail probability is determined by a small number of simulations 
and, therefore, has a relatively large error bar. In other words, 
the tails of the distributions of the different statistics are quite 
sparsely sampled. 

We will distinguish between the small (R < 10'), intermedi- 
ate (10' < R < 500') and the large (R > 500') scale regimes. Let 
us focus on the three statistics independently. We will highlight 
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Fig. 12. Difference of the normalized MFs obtained from the data with respect to the expected value of the null hypothesis for 
several sky coverages. The SMICA map is considered. From left to right and top to bottom: Area, Contour, Genus and Cluster- The 
grey bands represent the 1 and 2<x dispersions around zero, based on realistic Planck simulations including lensing, for / S k y = 0.23. 
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Fig. 13. Standard deviation (left), skewness (centre) and kurtosis (right) of the SMHW coefficients as a function of the wavelet 
scaled. Results are given for the four Planck CMB maps (green: Commander -Ruler, light-blue: NILC; red: SEVEM; yellow: SMICA). 
Modified upper tail probabilities (mUTP, see text for details) are obtained by comparing with 1000 simulations processed through 
the component separation pipelines. Squares represent modified upper tail probabilities that correspond to an actual upper tail 
probability above 0.5; diamonds represent upper tail probabilities below 0.5. 



the most important features and, afterwards, we will try to find 
an explanation for them: 



- On the smallest scales, the four CMB maps show a disper- 
sion in SMHW coefficients significantly larger than seen in 
the simulations. However on larger scales, the dispersion is 
systematically below the median of the simulations and, on 
scales of R « 5°, the modified upper tail probability is ap- 
proximately 0.015. 



- Regarding the skewness, all four maps yield a value that is 
significantly lower (with a modified upper tail probability 
of around 0.004) than expected from the simulations in the 
small scale regime (except for the smallest one, where the 
deviation is around 0.07). The rest of the scales are fairly 
compatible with the null hypothesis. 

- The kurtosis is also smaller than expected in the small scale 
regime. Overall, the modified upper tail probability is about 
0.03. At scales of around 300 r , an anomalously large value 
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(modified upper tail probability of approximately 0.01) is 
found. 

These re sults are compatible with the valu es reported for 
WMAP data flVielva et al.|2004l|Cruz et al.|20 05 ), over the scales 
common to both experiments (i.e., R > 10'). In particular, the 
large value of the kurtos is has been associated with the Cold 
Spot ( [Vielva et al. 2004| ). We wil return to this topic specifi- 
cally in Sect. |5.8| The low variance of th e wavelet coefficient s 
was previously seen in Viel va£taL] ( [2004] ) ; |Wiaux et al.| ( [2008] ) . 
In addition, the low dispersion at scales above a few degrees 
is likely to be related to the low variance anomaly detected in 
WMAP flMonteserm eTaTl|2008l ICruz et al.||20TT) , that is also 
seen in the Planck data (see Sect. 4.1 ). 



We have also studied the robustness of the results for differ- 
ent masking scenarios. In particular, we have investigated varia- 
tions in the results when we adopt, as auxiliary masks to define 
the exclusion masks, the two CG70 and CG60 masks removing 
30% and 40% of the sky, respectively. Note that the auxiliary 
masks obtained from the U73 mask already cut around 20% of 
the sky. The corresponding results for the SMICA map are pre- 
sented in Fig. 14 The conclusions are similar for the other CMB 
maps. For the dispersion of the wavelet coefficients, we do not 
notice any significant change in the anomalously high value ob- 
tained for the SMICA map at the smallest scales However, some 
changes are observed at larger scales. In this regime, it seems 
that the most significant deviation occurs for the CG70 mask 
(modified upper tail probability of around 0.005), whereas simi- 
lar (and slightly less significant) modified upper tail probabilities 
are obtained for both the U73 (modified upper tail probability of 
approximately0.015) and the CG60 (modified upper tail proba- 
bility of about 0.01) masks. A possible explanation for this be- 
haviour would be that a less restrictive mask admits some resid- 
ual contamination from Galactic foregrounds, thus increasing 
the dispersion of the wavelet coefficients, and artificially increas- 
ing their inconsistency with the null hypothesis. In principle, the 
larger the Galactic cut, the lower would be the dispersion of the 
wavelet coefficients (assuming that some residual contamination 
of the Galactic foregrounds is left) and, therefore, the smaller the 
upper tail probability. However, as we already said, the modified 
upper tail probability for the dispersion is higher for the CG60 
mask than for the CG70 mask. This apparent contradiction could 
be resolved by accounting for the larger sampling variance for 
smaller areas, that would result in a lower significance for the 
anomaly. 

The anomalous kurtosis at scales of R « 300' shows an over- 
all stable modified upper tail probability of around 0.01 — 0.03. 
In the small scale regime, the differences are better defined: the 
smaller the mask, the more significant the deviation (character- 
ized by the low value of the kurtosis). In particular, the mod- 
ified upper tail probability associated with the CG60 mask is 
0.001, around 0.009 for CG70, and approximately 0.03 for the 
U73 mask. A similar pattern is also observed for the skewness 
on these scales, although the three masks results in more similar 
upper tail probabilities, between around 0.001 and 0.007 (except 
for the smallest scale). 

It is therefore clear that there is some inconsistency between 
the CMB data and the corresponding simulations. On interme- 
diate scales, both the low dispersion and the high kurtosis could 
be related to previously known anomalies: the low variance and 
the Cold Spot. On the smallest scales, the three statistics report 
a low upper tail probability independently of the mask coverage 
— it is important to determine what this inconsistency is due to. 
Besides the possibility that it is an intrinsic cosmological sig- 



nal, the non-Gaussianity could be caused either by instrumental 
systematics or residual foreground contamination. 

In the former case, we have considered whether the origin 
of the signal could be related to properties of the noise that are 
inadequately modelled by the simulations. In particular, we have 
studied the statistical properties of the half-ring half-difference 
maps generated by the four component separation algorithms as 
proxies for the the noise present in the CMB maps. Although 
in detail there are some discrepancies between these noise es- 
timates and the simulated ones, they are not compatible with 
the inconsistencies observed between the CMB map and sim- 
ulations. Therefore, a systematic effect associated with the in- 
strumental noise does not provide a satisfactory explanation for 
the small-scale deviations. 

In the latter case, an obvious candidate is due to the con- 
tribution from residual unresolved point sources in the clean 
CMB maps. Although the brightest point sources are masked, 
and the component separation process itself can suppress the am- 
plitude of the unresolved b ackground of point sources, some si g- 
nal will remain. Indeed, in Planck Collaborati on XXlV| ( |2013| ) it 
has been determined that the bispectrum of this contribution is 
clearly detected in the four CMB Planck maps, at a significance 
in excess of 4<x. In addition, the dispersion of the wavelet coeffi- 
cients is higher than expected, which is also compatible with the 
presence of an additional signal. We therefore consider this as 
the most likely non-CMB explanation for the observed signal. 

4.6. Bispectrum 

The CMB bispectrum is the three point correlator of the at m co- 
efficients, 



(24) 



In this paper, we focus on the bispectrum reconstruction as 
a blind test of non-Gaussianity. Therefore, we assume we are 
seeking a non-trivial bispectrum that has arisen through a physi- 
cal process which is statistically isotropic, that is, we can employ 
the angle-averaged bispectrum B^^, 



m,2m 3 ' 



- / J n £ ] { 2 t 3 &m ] m 2 m 3 B m, 

where h^ 2 £ 3 is a geometrical factor, 



(2A + 1)(2£ 2 + 1)(2£ 3 + 1) / h h h 





An 



(25) 



(26) 



and Qmlmlml is me Gaunt integral, 



z7 m,\ mim 



= JdQ Y hm (n) Y hm2 (n) Y^ m3 (n) 

£\ £2 £3 
mi ni2 



= h^ 2i 



(27) 



with the usual Wigner-3 j symbol (m\m 2 m 3 \ It is more convenient 
to eliminate the geometrical factors entirely and to work with the 
reduced bispectrum which is defined as 



*w 3 - h hht B hhh 



(28) 



Note that the CMB bispectrum b^ 2 ^ 3 is defined on a tetrahedral 
domain of multipole triples {A £2^3} satisfying both a triangle 
condition and a limit given by the maximum resolution ^max 
of the experiment. A much more extensive introduction to the 
bispectrum can be found in Planck Collaboration XXIV ( 2013j ). 
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Fig. 14. Standard deviation (left), skewness (centre) and kurtosis (right) of the SMHW coefficients as a function of the wavelet scale 
R. Results are given for the SMICA CMB map. Several masking scenarios are compared: red: CG60 mask (cutting out 40% of the 
sky); green: CG70 mask (cutting out 30% of the sky); blue: U73 mask. The modified upper tail probabilities (mUTP) are defined in 
the text. 



Modal, wavelet and binned bispectrum estimators filter the 
CMB map with separable basis functions 



QijkVuh, = qi(t\)qj(h) qicWs) + perms , 



(29) 



to find the corresponding modal coefficients p^ (or p n because it 
is convenient to order the ijk with label n). For appropriately or- 
thonormalised basis functions Qijkifi, £2, £3), these coefficients 
can be used to reconstruct the CMB bispectrum through the 
signal-to-noise weighted expansion 



?W 3 



(30) 



This reconstruction method has been extensively validated, 
showing the accurate recovery of CMB bispectra from non- 
Gaussian simulated maps, and it has been applied to the WMAP 
seven year data to reco nstruct the full 3D CMB bispectrum 
(Fergu s son et al. 2010b). To quantify whether or not there is 
a model-independent deviation from Gaussianity, we can con- 
sider the total integrated bispectrum. By summing over all mul- 
tipoles, we can define an integrated nonlinearity parameter Fnl 
which, with the orthonormal modal decomposition (30), be- 
comes (Fergu sson et al.|2010b| ) 



NL 



1 h 2 b 2 

= A/2 Zj' 



loc £■ 



Tjijkfiijk 2 



(31) 



where A^ oc is the normalisation for the local /nl - 1 model (with 
coefficients ojjjjp. For ideal Gaussian CMB maps, the quantity 
should obey a x 2 -distribution with a mean given by the 
number of degrees of freedom (the modes) \i - n max and with 
a variance cr 2 = 2n max . Assuming that the three-point correla- 
tor is the leading non-Gaussian contribution, then Fnl provides 
a blind test for the presence of any integrated CMB bispectrum 
(once the expected two-point term is subtracted). We note that 
his is less sensitive than targeted searches for particular bispec- 
trum shapes. 

First, we discuss reconstructions from the modal estimator 
which has passed successfully through the full suite of non- 
Gaussian bispectrum validation tests (for further details about 
bispectrum estimators, see |Planck Co llaboration XXIV [2013 



We have applied this to the Planck temperature maps for the 
foreground- separation techniques NILC SEVEM and SMICA, us- 
ing two alternative sets of hybrid basis functions in order to 
cross-check results and identify particular signals. These are 




100 150 200 
Mode number n 



300 



Fig. 15. Planck recovered bispectrum coefficients p„ for the 
mode expansion ( [30} using hybrid Fourier modes (augmented 
with local and ISW modes). There is remarkable consistency be- 
tween results from the different component separation methods, 
NILC, SEVEM, and SMICA. The variance from simulated noise 
maps is nearly constant for each of the 300 modes, with the av- 
erage ±1<t variation shown in red. 



Fourier modes (n m2LX = 300) augmented a local SW mode and the 
separable ISW modes and a hybrid po lynomial/local basis with 
ftmax = 600, previously described in Ferg usson et al.| (f2010a). 
These basis function sets ensured excellent correlation with pri- 
mordial modes and the ISW signal. As with all the other bispec- 
trum analyses based on spherical harmonic coefficients, we used 
the U73 mask to which we applied inpainting. Together with 
the foreground separated maps, noise simulations were provided 
which were used to calibrate the estimator's linear correction 
term and to determine the variance. 

The modal coefficients p„ extracted from the Planck NILC, 
SEVEM, and SMICA maps are shown in Fig. [15] for the hybrid 
Fourier basis. These amplitudes show remarkable consistency 
between the different maps, with shape cross-correlations bet- 
ter than 96% and the overall amplitudes to within 7% agree- 
ment. This demonstrates that the indendent foreground separa- 
tion techniques do not appear to be introducing spurious non- 
Gaussianity. The p„ coefficients have a roughly constant vari- 
ance, so anomalously large modes can be easily identified. For 
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Fig. 16. Cumulative sum of orthonormal mode contributions J3„ 
to the total integrated bispectrum F[ 
ative quantity F 2 L = F 2 L - 1 NL M „^ , NL 

mean obtained from 200 CMB Gaussian maps and the standard 
deviation is the red line. A hybrid polynomial basis n max = 600 
is employed in the signal-dominated region £ < 1500. This^ 2 - 
test for the independent modes is cumulatively consistent with 
Gaussianity. 



example, we have subtracted the expected ISW signal and the es- 
timated point source contributions, explaining the large signal at 



low n. The corresponding quantity F 2 defined in (|31|), that can 



be seen in Fig. 16 shows consistency with the nulTTiypothesis. 
Using the modal expansion ( [30] ), we have reconstructed the full 
3D Planck bispectrum which is illustrated in Fig. [XT] for SMICA 



(large) but also NILC and SEVEM; the reconstructions are visu- 
ally indistinguishable. There are some striking features evident, 
notably the presence of a significant ISW modal contribution in 
the squeezed limit along the edges of the tetrapyd which has an 
oscillatory and flattened appearance. At large multipoles £ ap- 
proaching £ max = 2000, there is increased randomness in the 
reconstruction due to the rise in experimental noise and some 
evidence for a residual point source contribution. For the present 
Planck estimator configurations, the modal bispectrum estimator 
is more democratic, that is, it is capable of resolving the large-^ 
contributions near £ m?LX seen in Fig. [15] and not only the multi- 
poles associated with primordial models. 
In Fig. 
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we show a comparison of the £ < 500 Planck 
bispectrum signal and that reconstructed from the WMAP seven- 
year data ( Fergus so n et al.|2010b| ). Here for consistency we show 
the Planck signal from the second polynomial basis, since poly- 
nomials were used in the original WMAP1 analysis. The Planck 
signal pattern correlates well with the WMAP bispectrum ob- 
tained previously, despite the different domains used for the 
modal analysis of the two different experiments. 

Similarly to the modal bispectrum, a wavelet decompo- 
sition can be used to reconstruct the bispectrum. Here we 
use the continuous , non-orthogonal Spherical M exican Hat 
Wavelet (SMHW, |Martmez-Gonzalez et aT1|2002| ). Cubic mo- 
ments qijk are defined in terms of the SMHW coefficients for 
three different angular scales Ri, Rj, Rk (Curto et al. 2009b|a 
|20T0l|2011a|bl 



qtjk 



= -— f 

47T CTiCTjCTk J 



dnw(Ri, ri)w(Rj, ri)w(Rk, n) 



(32) 



where <x; is the dispersion of the wavelet coefficient map 
w(Ri, n). Considering the covariance matrix of the q^k moments, 



Table 16. x 2 statistics based on the wavelet bispectrum recon- 
struction yi statistics for the foreground cleaneddatadata map. 
Considered data map: combined map cleaned with C-R, NILC, 
SEVEM, and SMICA. 



Method 


/^data 


DOF 


<x 2 ) 


(T 




C-R 


874 


690 


740 


87 


0.074 


NILC 


883 


682 


731 


83 


0.045 


SEVEM 


858 


682 


731 


83 


0.070 


SMICA 


. . 878 


682 


732 


83 


0.058 



C = <qq T ), and its eigenvector decomposition, C = RDR T , with 
R the eigenvector matrix and D the eigenvalue matrix, a new set 
of quantities y = D 1/2 R T q is defined. Considering the decorrela- 
tion produced by the convolution of the SMHW on the temper- 
ature anisotropies and applying the central limit theorem to the 
averages defined in Eq. [32] then the quantities are expected 
to have a nearly Gaussian distribution. Therefore, the y quanti- 
ties are nearly mu ltinormal and satisfy (yy T ) = I and (y) = 
dCurto et al.|20TTa]) . 

We have computed this reconstruction using the Planck data 
and compared with the null hypothesis (Gaussian Planck sim- 
ulations). The considered data map is the resulting map af- 
ter foreground cleaning based on different cleaning procedures: 
Commander -Ruler, NILC, SEVEM, and SMICA. The mask used is 
the U73 one (contrary to the modal reconstruction, no inpainting 
of the point sources is made in this case). In Fig. [19] the y statis- 
tics corresponding to the Planck data are plotted and compared 
with the 3cr error-bars obtained with Planck Gaussian simula- 
tions. From the list of different q ^ statistics corresponding to 
the 16 angular scales described in Planck Collaboration XXIV 
(2013 ), there are 11, 4, 3, 3 statistics with \y t \ > 3 (corresponding 
to Commander-Ruler, NILC, SEVEM, and SMICA respectively). 
The error-bars are obtained with Planck simulations for each 
type of component separation cleaned map. The error-bars of the 
yt statistics for low indices i are associated to large scales where 
the q statistics have a less Gaussian-like shape. The y statistics 



are combined into a x t est a fter a principal com ponent analy- 



sis with a threshold of 10 12 ( Curto et aT]|20 1 1 a| ) and compared 
with the^ 2 statistics obtained from Planck Gaussian simulations 
for each type of component separation method (see Table 16). 
The^ 2 statistic corresponding to the data is compatible with the 
values obtained from Gaussian simulations according to the cu 



16 
3e- 



mulative probability P(x 2 > Xd ats )i as can be seen in Table 
Therefore the wavelet bispectrum reconstruction does not 
tect a significant amplitude of bispectrum in the considered data 
maps. Details on the constraints on the amplitude of different 
bispectrum shapes are presented in |Planck "Co llaboration XXIV 



5. Intriguing inconsistencies - WMAP anomalies 
revisited. 

In the previous section, we have established that the Planck data 
shows little evidence for non-Gaussianity beyond that expected 
due to the ISW-lensing effect (which is accounted for directly by 
simulations), and contributions from residual unresolved point 
sources. The exceptions are on large-angular scales where fea- 
tures consistent with various anomalies previously seen in the 
WMAP data have been observed. In this section, we explic- 
itly consider several of the most important anomalies detected 
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500 2000 



Fig. 17. Full 3D CMB bispectrum recovered from the Planck foreground-separated maps (modes illustrated in fig. 



15 



including 

SMICA (left), NILC (centre) and SEVEM (right). These are plotted in three-dimensions with multipole coordinates [£\, £2, £3}; the 
triangle condition restricts the bispectrum to a tetrahedral domain out to the experimental resolution limit I < £ max = 2000. Several 
density contours are plotted with red positive and blue negative. The bispectra from different component- separation methods are 
almost indistinguishable with the same features also appearing in Fourier and polynomial expansions. Note the central and flattened 
features for i < 1200 and also the oscillating CMB ISW lensing signal in the squeezed limit along the edges of the tetrapyd. 




Fig. 1 8. Comparison between t he WMAP seven-year bispectrum signal 
(left) ([Fergusson et al. (2010b )) and the \ow-£ signal of Planck (right) 
reconstructed from the SMICA foreground-separated map (in both cases 
using polynomial modes). The same basic patterns are observed in both 
bispectra, including an apparent central 'oscillatory' feature. 
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in the WMAP data, namely the quadrupole-octopole alignment 
(Sect. |5.1|), th e low variance (Sect. |5.2| ) , hem ispherical asymme- 
try (Sect. 5.3 }, pha se correlations (Sect. |5.4[ ), dipolar power mod- 
ulation (Sect. [53}, generalized power modulation (Sect. |5.6| ), 
parity asymmetry (Sect. |5.7( ) and the Cold Spot (Sect. |5.8( ). Each 
of these anomalies may represent different violations of the fun- 
damental properties of isotropy and/or Gaussianity of the CMB 
data which are assumed in the estimation of the CMB power 
spectrum. 

There is an ongoing debate about the significance of these 
anomalies in the literature. A critical issue relates to the role of 
a posteriori choices — whether interesting features in the data 
bias the choice of statistical test or if arbitrary choices in the sub- 
sequent data analysis enhance the significance of the features. 
Indeed, the WMAP team ( [Bennett et al.|20 11 1 ) contends that the 
anomalies are significantly over-interpreted due to such selec- 
tions, whilst other authors claim highly significant and robust 
detections. Therefore, care must be taken to address the issue, 
since our analyses are necessarily follow up tests of the previ- 
ous WMAP investigations. However a careful and fair statistical 
treatment can allow us to study possible links among the anoma- 
lies and to search for a physical interpretation. 



Fig. 19. The wavelet bispectrum reconstruction y t statistics for 
the foreground cleaned Planck data map. Considered data map: 
combined map cleaned with C-R, NILC, SEVEM and SMICA. The 
solid yellow lines represent the 3<x error-bars for SMICA (similar 
error-bars are obtained for C-R, NILC, and SEVEM maps). 

5.1. Mode alignment 

Tegmark et al.| ( |2003| ) first reported a significant alignment be- 
tween the orientation of the quadrupole and the octopole in the 
WMAP first year temperature data. We study this quadrupole- 
octopole alignment in the Planck data using the max imization 
of the angular momentum dispersion as described in de Oliveira- 
|Costa et aT| ([2004). Specifically, we determine the orientation of 
the multipoles by finding the axis n around which the angular 
momentum dispersion 



1 



m 2 \a £m (n)\ 2 



(33) 



is maximized. Here, ai m (ri) denotes the spherical harmonics co- 
efficients of the CMB map in a rotated coordinate system with 
its z-axis in the fi-direction. This definition of the multipole- 
orientation has been devised for planar multipoles and is sim- 
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ply the direction perpendicular to the plane in which most of the 
power of the multipole lies. It is thus intuitive and easy to use. 
Note that the value of the statistic in Eq. ( [33] ) is the same for 
-n as for #i, i.e. the multipole orientation is defined only up to a 
sign. 

An al ternative method, based on the multipo l e vector decom- 
positi on ([Copi et al.|2004[|Schwarz et al.|2004| |Bielewicz et al. 
2005 ; Bielewicz & Riazuelo 2009) of the data has also been used 
to verify the robustness of the results presented here, and excel- 
lent consistency is found. 

Residual foregrounds (mostly on the Galactic plane) present 
in the four Planck CMB estimates could influence the recon- 
struction of the low-order multipoles. However, when a mask 
applied, the resulting mode-coupling can also affect the recon- 
struction of the low-£ multipoles. We therefore utilise Wiener 
filtered maps computed from the data to which the U73 mask is 
applied. Specifically, we utilise the same implementation of the 
Wiener filter as used in Planck Collaboration P09A ( 20 13]) i.e., a 
messenger method as described by Eisner & Wan delt| ([201 3|). 
A direct inversion method fo r masked data (Efstathiou 2004; 
Biel ewicz et al. 2004 2013| ) is a possible alternative, but the 



Wiener filtered maps result in a significantly smaller uncertainty 
in the reconstructed orientation of the multipoles. 

We then search for the preferred orientation by explicitly ro- 
tating the CMB map such that the z-axis pierces the centre of all 
the low resolution pixels defined at N S ^ Q = 16, and then subse- 
quently refine the search by using an A^e = 2048 map. The an- 
gular resolution for the orientation of the multipoles is thus given 
by the distance between the pixel centers of the N S ^ Q = 2048 
map, which is of order 1.94'. Figure 20 shows the Wiener fil- 
tered SMICA CMB sky, with the corresponding reconstruction of 
the quadrupole and octopole moments. The reconstructed orien- 
tations are quite robust with respect to the component separation 
method used for reconstructing the CMB. The significance of 
the alignment between the quadrupole and the octopole is as- 
sessed from the scalar product of their orientations, compared to 
values derived from the standard set of 1000 MC simulations. 
The orientation, the angular distance the scalar-product between 
quadrupole and octopole, and the probability of at least such an 
alignment to occur in an isotropic universe are summarised in 
Table[T7]for each CMB map. 

We find that, depending on the component separation 
method, the quadrupole and octopole orientations are mis- 
aligned by an amount betwe en 9° and 13°. This is larger than 
the 3° reported recently by Be nnett et al.| d2012) for the 9- 
year WMAP ILC map. In consequence, our significance of the 
quadrupole-octopole alignment is substantially smaller than for 
the WMAP data, falling to almost 98% confidence level for 
the Commander -Ruler and SEVEM maps and 96.7% confidence 
level for the NILC map. 

5.2. Variance, skewness and kurtosis anomalies 

A low value for the variance on the CMB sky was previously ob- 
served in th e WMAP data by |Monteserin et a l . ( 20 08 ) and |Cruz 
|et al.| ( [20TT] ), and confirmed for Planck in Sect. |4.1| Furtherm ore, 
the effect has also been seen in the wavelet analysis of Sect. 4^5_ 
where the variance of the SMHW coefficients is low at scales 



between 400 and 600 arcmin (Fig. [13]). In addition, anomalous 
behaviour was also observed for the skewness and kurtosis in 
low resolution maps at N s ^ e = 16. Here, we reassess these re- 
sults and determine their robustness to masking and data selec- 
tion. The former will allow us to determine whether a particular 





Fig. 20. Upper: The Wiener filtered SMICA CMB sky (temper- 
ature range ± 400 jjK). Middle: the derived quadrupole (tem- 
perature range ±35 yuK). Lower: the derived octopole (temper- 
ature range ± 35 yuK). Cross and star signs indicate axes of the 
quadrupole and octopole, respectively, around which the angular 
momentum dispersion is maximized. 



region is causing the anomalous behaviour, whilst the latter can 
establish whether foreground residuals could be responsible. 



Table 18 and Fig. 21 present the results for the variance, 
skewness and kurtosis determined from the four CMB maps 
with the U73, CL58 and CL37 masks applied. Results are also 
computed for data within the ecliptic hemispheres surviving 
the U73 mask. The variance is low in all cases, with only 
small differences in significance observed for the different maps. 
Interestingly, the low variance seems to be localised in the north- 
ern ecliptic hemisphere. Conversely, anomalous values for the 
skewness and kurtosis are only apparent for the southern ecliptic 
hemisphere. 

Since these results might be indicative of the presence of 
Galactic foreground residuals near the Galactic plane, we anal- 
yse the frequency dependence of the statistics as summarised in 
Table 19 and Fig. 22 The variance shows little frequency depen- 
dence for the considered masks and regions, whereas the skew- 
ness and kurtosis show a moderate frequency dependence when 
the U73 mask is applied, as also seen for the N-pdf in Sect. [42 



20 



Planck Collaboration: Isotropy and statistics 



Table 17. Orientation of the low multipoles extracted from the different component separated CMB maps, obtained from maximizing 
the angular momentum dispersion. The second last column gives the absolute value of the scalar-product between the orientation 
vectors of the quadrupole and the octopole. In an isotropic universe, the latter is uniformly distributed on the interval [0,1]. The last 
column gives the probability of such an alignment (or stronger than that) to occur. 



Method (l,b) quadrupole [°] (l,b) octopole [°] ang. distance [°] scalar-product probability 



C-R (228.2,60.3) (246.1,66.0) 9.80 0.985 0.019 

NILC (241.3,77.3) (241.7,64.2) 13.1 0.974 0.033 

SEVEM (242.4,73.8) (245.6,64.8) 9.08 0.988 0.016 

SMICA (238.5,76.6) (239.0,64.3) 12.3 0.977 0.032 



|Cruz et aL] (2011) found that a small region of the sky localised 
to both the ecliptic and Galactic south and near to the Galactic 
plane (their so-called gplO region) exhibited particularly high 
variance. Thus, since the skewness is negative, we consider a 
prominent cold spot at (b = -8°, / = 32°), partially masked by 
the Galactic plane. However, when masking the seven coldest 
pixels of the spot, the significance of the skewness and kurto- 
sis drops only slightly, with lower tail probabilities of approx- 
imately 0.03 and 0.93 respectively. If the whole gplO region 
(/sky = 7%) is masked, the skewness and kurtosis drop dras- 
tically and have lower tail probabilities of approximately 0.30 
and 0.50 respectively, whereas the variance is highly significant 
since none of the 1000 simulations has a variance below the data. 
In order to check the possible leakage of Galactic contamina- 
tion due to the Gaussian smoothing applied to the low resolution 
data, we rep eated our calculations for the Wiener filtered maps 
used in Sect. |5.1| but found little variation to the existing results. 
Therefore, it is unlikely that any leakage impacts the estimators. 

The incompatibility of the observed variance with simula- 
tions based on a cosmological model that has been determined 
from the same data set might appear puzzling at first, but can 
be understood as follows. The map-based variance is dominated 
by contributions from large angular scales on the sky, whilst 
the cosmological parameter fits are relatively insensitive to these 
low-order ^-modes, and are instead largely dominated by scales 
corresponding to £ > 50. Thus, the best-fit spectrum in the con- 
text of a 6-parameter ACDM model can have a mismatch with 
the data on these scales, so that the corresponding simulations 
will not adequately capture the dearth of power at low-£. The re- 
sults presented here do indeed imply that the large-angular scale 
power is low relative to the fiducial sky model. In fact, when 
subtracting the quadrupole and octupole from both the data and 
simulations outside the U73 mask, the results are more consis- 
tent. In this case, the lower tail probabilities for the variance, 
skewness and kurtosis are 0.192, 0.637 and . 792 re spectively. 
This result was already found in Cruz et al. ( 2011| ). It is then 
plausible that the low multipole alignment could have the same 
cause as the anomalies considered here. However, when sub- 
tracting the quadrupole and octupole outside the CL58 mask, 
the lower tail probability for the low variance is 0.036, which 
remains rather low. The connection with the very low power in 
the ecliptic northern hemisphere also remains to be explored. 



5.3. Hemispherical Asymmetry 



Table 18. Lower tail probablity for the variance, skewness and 
kurtosis at A^e = 16, using different masks. 



In |Eriksen et aL] (2004a) and |Hansen et aT] (2004) it was dis- 
covered that the angular power spectrum of the first year WMAP 
data, when estimated locally at different positions on the sphere, 
appears not to be isotropic. In particular, the power spectrum 
calculated for a hemisphere centered at (0,0) = (110°, 237°) 
(in Galactic co-latitude and longitude) was larger than when 



Mask 



C-R NILC SEVEM SMICA 



Variance 



U73,/ sky =78% 0.019 0.017 0.014 0.019 

CL58,/ sky =58% 0.004 0.003 0.003 0.003 

CL37, / sky = 37% 0.028 0.017 0.018 0.016 

Ecliptic North, / sky =39% 0.001 0.001 0.001 0.002 

Ecliptic South, / sky =39% 0.464 0.479 0.454 0.490 



Skewness 



U73,/ sky =78% 0.016 0.015 0.023 0.012 

CL58, / sky = 58% 0.208 0.139 0.162 0.147 

CL37, / sky = 37% 0.517 0.467 0.503 0.469 

Ecliptic North, / sky = 39% 0.502 0.526 0.526 0.521 

Ecliptic South, / sky =39% 0.004 0.006 0.008 0.004 



Kurtosis 



U73,/ sky =78% 0.972 0.973 0.966 0.982 

CL58,/ Sky =58% 0.630 0.726 0.711 0.711 

CL37, / sky = 37% 0.069 0.135 0.130 0.124 

Ecliptic North, / sky = 39% 0.094 0.229 0.196 0.245 

Ecliptic South, / sky = 39% 0.933 0.916 0.886 0.948 



calculated in the opposite hemisphere over the multipole range 
£ = 2-40. Simultaneously, Park (2004 ) also presented evidence 
for the existence of such hemispherical asymmetry — in which a 
particular statistical measure is considered to change discontin- 
uously between two hemispheres on the sky — with the appli- 
cation of Minkowski functionals to the WMAP data. Since the 
preferred direction of Eriksen et al.| (2004a| ) lies close to the 
ecliptic plane, it was also demonstrated that the large-angular 
scale N-point correlation functions showed a difference in be- 
haviour when computed on ecliptic hemispheres. Many studies 
have subsequently been undertaken fo cusing on hemispher es in 
the ecliptic coordinate system, with |Schwarz et al.| (2004) par- 
ticularly emphasizing the connection. Hemispherical asymme- 
try has also been seen with other measures o f non-Gaussianity 
(Eriksen et al.|2004c| [20051 |Rath et al.|2007a| ). 

Here we repeat the analysis of Eri ksen et al.| (2005) on the 
Planck component separated data, smoothed and then down- 
graded to_/V S ide = 64 as described in Sect. [2] As already noted 
in Sect. |4.3| the results for the low resolution maps are the most 
deviant relative to the MC simulations based on the Planck fidu- 
cial model. 

The N-point correlation functions computed on the northern 
and southern hemispheres determined in the ecliptic coordinate 
frame and using the U73 mask are shown in Fig. 23 The cor- 
relation functions for the four Planck maps are very consistent, 
and the observed behaviour is in agreement with that seen in 
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Fig. 21. Variance, skewness and kurtosis at N^q = 16, for the 
U73 mask, CL58, CL37, ecliptic North and ecliptic South (from 
top to bottom). The different lines represent the four component 
separation methods C-R (green), NILC (blue), SEVEM (red), and 
SMICA (orange). 



Fig. 22. Variance, skewness and kurtosis at A^e = 16, for the 
U73 mask, CL58, CL37, ecliptic North and ecliptic South (from 
top to bottom). The different lines represent the four considered 
frequencies, namely 70 GHz (green), 100 GHz (blue), 143 GHz 
(red), and 217 GHz (orange). 



Table 19. Frequency dependence of the lower tail probablity for 
the variance skewness and kurtosis at A^de = 16, using different 
masks. 



Mask 



70 GHz 100 GHz 143 GHz 217 GHz 



Variance 



U73,/ sky =78% 0.019 0.013 0.014 0.016 

CL58,/ sky =58% 0.003 0.003 0.003 0.003 

CL37, / sky = 37% 0.024 0.020 0.018 0.020 

Ecliptic North, / sky = 39% . 0.001 0.002 0.001 0.001 

Ecliptic South, / sky = 39% . 0.446 0.436 0.455 0.455 

Skewness 

U73,/ sky =78% 0.045 0.016 0.024 0.015 

CL58,/ sky =58% 0.254 0.205 0.162 0.157 

CL37, / sky = 37% 0.503 0.471 0.468 0.515 

Ecliptic North, / sky = 39% . 0.505 0.447 0.541 0.352 

Ecliptic South, / sky = 39% . 0.015 0.006 0.009 0.006 

Kurtosis 

U73, / sky =78% 0.962 0.981 0.965 0.974 

CL58, / sky = 58% 0.619 0.684 0.710 0.725 

CL37,/ sky =37% 0.114 0.091 0.130 0.121 

Ecliptic North, / sky = 39% . 0.180 0.096 0.203 0.180 

Ecliptic South, / sky = 39% . 0.902 0.920 0.882 0.909 



the WMAP data ( [Eriksen et aL] [2004a). Specifically, the north- 
ern hemisphere correlation functions are featureless (both the 
three- and four-point functions lie very close to zero), whereas 
the southern hemisphere functions exhibit a level of structure 
that is in good agreement with the confidence regions computed 
from the Gaussian simulations. 

The probabilities of obtaining a value for the x 2 statistic for 
the Planck fiducial ACDM model at least as large as the ob- 



served values are presented in Table [20] The probabilities for the 
3 -point and 4-point functions in the northern Ecliptic hemisphere 
are especially large, and in the case of the pseudo-collapsed con- 
figuration all simulations yielded a larger than observed value of 
the^ 2 . Nominally, this value is even more remarkable than found 
with the WMAP data (Eriks en et al.|2004a| , although to interpret 
it correctly one has to keep in mind that the analysis presented 
here is an example of a posteriori statistic. Specific choices have 
been made about the smoothing scale used for downgrading the 
data, and, in particular, for the reference direction that defines 
the hemispheres. This will tend to overestimate the significance 
of the results. Nevetheless, the observed properties of the Planck 
data are consistent with a remarkable lack of power in a direc- 
tion towards the north ecliptic pole, consistent with the simpler 



one-point statistics presented in Sect. 5.2 



5.4. Phase correlations 

Previous studies using the methods of scaling indices and surro- 
gates and based on the WMAP three-, five- and seven-year data 
(Rat h et al.|[2009] |20TT| |Rossmanith eTaLl|20T2l |Modest et al.l 
2013} showed significant evidence for intrinsic phase correla- 
tions at low € values in the CMB. The signal was demonstarted 
to be robust with respect to the WMAP data release, to the com- 
ponent separation methods and to the selected test statistics. In 
this section we apply these methods to the Planck component 
separated data sets. 

The scaling index method represents one way to estimate 
the local scaling properties of a point set in an arbitrary d- 
dimensional embedding space. The technique provides the pos- 
sibility to reveal local structural characteristics of a given point 
distribution. A number of analyses have used scaling indices 
to test the Gaussian nature and statistical isotropy of the CMB 
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Fig. 23. The 2-point (upper left), pseudo-collapsed (upper right), equilateral 3-point (lower left) and rhombic 4-point (lower right) 
correlation functions (Af s ide = 64). Correlation functions are shown for the analysis performed on northern (blue) and southern 
(red) hemispheres determined in the Ecliptic coordinate frame. The shaded dark and light grey bands indicate the 68% and 95% 
confidence regions, respectively. 



as represented by the WMAP data (|Rath et al.||2007a| [2009; 
IRossmanith et al.|2009a||Rath et aL|20TTF ~ 

In general, the method is a mapping that calculates, for each 
member p u i = 1, . . . , N p { x of a point set P, a single value that de- 
pends on the spatial position of pi relative to the group of other 
points in its neighborhood, in which the point under consider- 
ation is embedded. A three-dimensional point set P is gener- 
ated for two-dimensional spherical CMB-data by transforming 
the temperature values T(0i, (pi) of each pixel to a radial jitter 
around a sphere of radius R at the position of the pixel centre 
(0j, (pi). For obtaining scaling indices the local weighted cumula- 
tive point distribution which is defined as 



piPur) = ^s r (d(j>i,pj)) 



(34) 



7=1 



with r describing the scaling range, while s r and d denote a shap- 
ing function and a distance measure, respectively, is calculated 
first. The scaling index a(pi, r) is then defined as the logarithmic 
derivative of p(p i9 r) with respect to r: 



®(Pi, r) 



dlogp(pi,r) 
dlogr 



(35) 



Using a quadratic gaussian shaping function s r (x) = e~^ 2 and 
an isotropic euclidian norm d(pi, pj) = \\pi ~Pj\\ as distance mea- 
sure, one obtains the following analytic formula for the scaling 
indices 



a(Pi, r) = 



(36) 



where we use the abbreviation dy = d(pi,pj). As should be 
clear from equation (36), the calculation of scaling indices de- 
pends on the scale parameter r. Ten scaling range parameters 
r k = 0.05, 0.1,..., 0.5, k = 1,2,.. 



10 in the notation of Rath 
et al. (2007a]) are used in this analysis. In order to calculate scal- 
ing indices on large scales as in previous studies, we couple the 
r-jitter a to via a - 0.5r^. The mean (a(r^)> and the standard 
deviation cr a ( rk ) derived from the full sky and from a set of 768 
rotated hemispheres are used to test for non-Gaussianity and de- 
viations from statistical isotropy. 

In order to quantify the significance of the scaling index 
results, and focus the study on the phase properties of the ob- 
served CMB sky, we utilize the method of surrogate maps ( Rath 
et al.||2009] ). Such a technique offers the unique possibility to 
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Table 20. Probabilities of obtaining values of the x 2 statistic for 
the Planck fiducial model at least as large as the observed values 
of the statistic for the Planck maps with resolution parameter 
Af side = 64 estimated using the C-R, NILC, SEVEM and SMICA 
methods. 





C-R NILC SEVEM 


SMICA 


Two-point function 


Northern Ecliptic 
Southern Ecliptic 


0.935 0.924 0.927 
. . 0.633 0.599 0.639 


0.932 
0.592 


Pseudo-collapsed three-point function 


Northern Ecliptic 
Southern Ecliptic 


1.000 1.000 1.000 
. . 0.349 0.310 0.381 


1.000 
0.301 


Equilateral three-point function 


Northern Ecliptic 
Southern Ecliptic . . 


0.996 0.999 0.999 
. . 0.627 0.644 0.678 


0.999 
0.656 


Rhombic four-point function 


Northern Ecliptic 
Southern Ecliptic . . 


0.999 0.999 0.999 
. . 0.559 0.548 0.574 


0.999 
0.553 



test for scale-dependent non-Gaussianity and deviations from 
isotropy in a completely model-independent ("blind") way. This 
self-consistency of the surrogate approach suppresses the sensi- 
tivity of the null tests to the assumed fiducial power spectrum. 
This is particularly pertinent given the potential mismatch of 
the Planck data to the fiducial sp ectrum on large-angular scales 
( Planck Co llaboration XV 2013 ). The statistical properties of a 
Gaussian random field on the sphere can be fully described by its 
two-point correlation function (or power spectrum) and exhibit 
Fourier phases that are independent and identically distributed 
(i.i.d.) and follow a uniform distribution in the interval [-7r,7r]. 
Thus, demonstrating the existence of Fourier phase correlations 
in CMB maps could indicate the presence of non-Gaussianities 
in the primordial density fluctuations. The possible presence of 
phase correlations is tested using the method of surrogates. 

However, the Gaussianity of the temperature distribution 
and the randomness of the set of Fourier phases of the map to 
be studied are a necessary pre-requisite for the application of 
the surrogate-generating algorithm. Therefore the following pre- 
processing steps are applied to generate a zero order surrogate 
map. First, the maps are remapped on to a Gaussian distribution 
in a rank-ordered way. Then we ensure the randomness of the 
set of Fourier phases by making a rank-ordered remapping of 
the phases on to a set of uniformly distributed ones. These two 
preprocessing steps only have marginal influence on the maps. 
Now, the set of surrogates to be used for assessing the statis- 
tical properties of the data sets can be generated by shuffling 
the phases in the space of the spherical harmonics while exactly 
preserving the modulus of the ai m . Moreover, by introducing a 
two-step shuffling scheme for previously specified ^-ranges, a 
scale-dependent analysis is made possible. It is worth noticing 
that while in all surrogate maps the modulus of the ai m is ex- 
actly preserved, null tests involving a comparison to an assumed 
fiducial power spectrum only preserve the Q values, which are 
average values of the \az m \ when summed over m. Thus, the lin- 
ear properties of the surrogate maps are more tightly constrained, 
and specifically kept constant, than in tests involving simulated 
maps generated on the basis of the Qs. 

So-called first and second order surrogate maps are then 
obtained as follows. We initally generate a first order surrogate 



map, in which any phase correlations for the scales that are 
not of interest are randomized by shuffling its phases (p£ m for 
£ £ Al = [Cin, 4iaxL < m < £. In a second step, N (N = 1000 
throughout these investigations) realizations of second order 
surrogate maps are generated from the first order surrogate map, 
in which the remaining phases <f>£ m with £ e A£, < m < £ are 
shuffled, while the previously randomized phases for the other 
scales are preserved. The generation of surrogates is always 
performed using the maps with the highest resolution, i.e., 
A^side = 2048. Given the evidence for anomalies on the largest 
angular scales, and to ensure consistency with the previous 
WMAP analyses, we perform dedicated scale-dependent tests 
for the scales defined by A£ = [2, 20]. 

Since the methodology in its simplest form requires the or- 
thonormality of the set of basis functions Y^ m , we apply the 
method to the full sky foreground-cleaned maps as obtained after 
component separation with Commander-Ruler, NILC, SEVEM 
and SMICA. For the selected ^-interval A£ = [2, 20], the genera- 
tion of the first order surrogate map removes the phase signature 
of the small scale residuals in the data and can be interpreted 
as a form of inpainting procedure for small masked patches in 
the Galactic plane. The differences between the first and second 
order surrogates are quantified by the <x-normalized deviation S 



S(Y) 



' surrol 



(^surro2) 



(37) 



with, Y = (a(rk)),cr a ( rk ),x 2 - Here, x 2 represents either a diago- 
nal combination of the the mean (a(rk)) and standard deviation 
o- a (r k ) at a certain scale or for the full scale-independent x 2 
statistics 



X 2 = (M-(M)) T C-\M-(M)), 



(38) 



where the test statistics to be combined are comprised in the 
vector M and C is obtained by cross correlating the elements of 
M . With the mean and the standard deviation as input for M we 
obtain x\ a) and xl- a statistics with M T = ((a(ri)), . . . , (a(no))) 
and M T = (cr a(nh . . . , cr a{no) ) respectively. 



Fig. 
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shows the Six 2 ) values for the set of rotated hemi- 
spheres for the SMICA map. Each pixel of the full sky map 
with a HEALPix resolution of = 8 specifies one of the 
768 S-values for a rotated hemisphere, where the pixel posi- 
tion indicates the orientation of the z-axis of the rotated coor- 
dinate system. We detect pronounced signatures for both non- 
Gaussianities and anisotropies. The results are consistent for 
Commander-Ruler, NILC, SEVEM and SMICA. 



In Fig. [25] the deviations S (Y) are displayed for the mean and 
standard deviation. We only show the results for the SMICA map. 
The other three maps yield very similar results. For all four maps 

the values for S ((a)) extend far beyond 3 for r = 0.2 0.25 

when rotated hemispheres are considered separately. Since the 
effect in the separate hemispheres goes in opposing directions, 
no signal is observed for the full sky. The results for the scale- 
independent^ 2 statistics are summarized in Table 21 The results 
suggest a highly significant detection of both non-Gaussianities 
and anisotropies in the Planck data, consistent with t hose ob- 
tained previously with WMAP data (for comparison see Modest 
|et al.|2013| ). 

We have also investigated whether the significance of the 
results depends on the choice of £ m \ n and £ m2LX . In particular, 
we have extended the range of interest to £ max = 30, and then 
considered three sub-intervals, A£ = [2, 10], A£ = [11,20] and 
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Fig. 24. Deviations Six 2 ) °f tne rotated hemispheres derived 
from a combination of the mean and the standard deviation of 
the scaling indices for the scale r$ determined from the SMICA 
map. 




Fig. 25. Deviations \S(r)\ for the SMICA map as a function of 
the scale parameter r for the full sky (black) and upper (red) 
and lower (blue) rotated hemispheres. The plus signs denote the 
results for the mean <ar(r^)>, the star-signs represent the standard 
deviation cr a ( rk y The dashed (dotted) line indicates the 1 (3) <x 
significance interval. 



Table 21. Deviations S and empirical probabilities p for the 
scale-independent ^-statistics derived from the C-R, NILC, 



SEVEM and SMICA maps. 



Full Sky Upper Lower 

Hemisphere Hemisphere 

(S/%) os/%) os/%) 



C-R,*< a> 0.86/82.6 4.21/99.7 3.18/99.0 

C-R,xl a 0.88/85.2 3.94/99.5 3.10/99.2 

NILC,^ 0.86/81.8 3.74/99.6 4.41 / >99.9 

NILC,^ 0.79/78.8 3.69/99.6 4.49/>99.9 

SEVEM, x\ a ) 0.00/58.0. 3.22/99.3 5.02/>99.9 

SEVEM, xl- a 0.05/60.8 3.20/99.0 5.11/99.9 

SMICA, X \ a) 0.75/80.1 3.80/99.8 4.70/99.8 

SMICA, xl- 0.01/54.4 3.64/99.3 4.81/>99.9 



At = [21, 30]. Over the full range, similar results are found to 
those from A£ = [2, 20], but at lower significance. This suggests 
that the inclusion of the phases in the interval At = [21, 30] di- 
lutes the signal, because there are no phase correlations in this 
Grange. This is corroborated by the fact that the first and sec- 



ond order surrogates generated specifically for this sub-interval 
cannot be distinguished. The results for the interval At = [2, 10] 
are quite consistent with those over A£ = [2,20], whereas for 
At = [11, 20] we find that the signature in the northern ecliptic 
hemisphere nearly vanishes. Conversely, in the southern eclip- 
tic hemisphere, on the other hand, the S -signal persists - espe- 
cially for regions covering the Cold Spot. It thus appears that 
the lowest /'-range is predominantly responsible for the detected 
hemispherical asymmetries detected in the spectrum of scaling 
indices, whilst the interemediate interval considered may have 
an association with the Cold Spot. It is certainly the case that 
scale-dependence is seen in the nature of the phase correlations 
present in the data. 

Since both, the modulus of the a^ m s for all £s and the phases 
(f)£ m for £ g A£ are exactly the same in the first and second or- 
der surrogates, one must infer that the pattern of hemispherical 
asymmetry in the S-maps can solely be attributed to phase corre- 
lations in the interval At. Thus, the analysis involving surrogate 
maps reveals that there are phase correlations at low I. 

5.5. Dipolar asymmetry 

In previous sections, we have seen evidence for a break in 
isotropy related to the discontinuous distribution of power in 
hemispheres on the sky. Bennett et al. ( 2011) distinguishes be- 
tween such an asymmetry and one where the CMB signal is 
modulated across the sky by a dipolar term. Studies of such a 
dipolar asymme try have been motivated by the phenomenolog- 
ical proposal of |Gordon et al.| ( [2005| ) that the power asymmetry 
could be described in terms of a multiplicative dipole modula- 
tion model. In addition, relativistic Doppler boosting due to our 
motion with respect to the CMB rest frame is expected to in- 
duce a dipolar modulation aligned with the CMB dipole at the 
<9(10~ 3 ) level; a statistically significant detection of this effect 
by Planck is presented by |Planck Collaboration XXVII| ( [20T3 ). 



5.5.1 . Power asymmetry 



In their analysis of the 5-year WMAP data, (Ha nsen et al.|2009| ) 
specifically searched for the dipolar modulation of power on 
the sky, In particular, a simple test was performed in which the 
power spectrum on discs was computed and binned into inde- 
pendent blocks of 100 multipoles from I - 2 to I - 600, then 
each block fitted for a dipolar asymmetry in the power distribu- 
tion. The 6 ^-ranges considered showed evidence of a consis- 
tent dipole direction, yet, from a set of 10000 simulations, none 
showed a similarly strong asymmetry. A further extension of the 
analysis introduced a model selection procedure taking into ac- 
count the statistical penalty for introducing an asymmetric model 
with additional parameters (direction of asymmetry, amplitude 
of asymmetry and asymmetric multipole ranges). Even in this 
case, the asymmetry was found to be highly significant for the 
whole range £ = 2- 600. 

However, such a procedure is highly expensive in terms of 
CPU-time. Given the higher sensitivity and angular resolution 
of the Planck data, we have therefore elected to focus on the 
simpler disc-based test, thus allowing us to probe further into a 
previously unexplored Grange. This should at least in part an- 
swer any a posteriori criticisms of the study. Since the analysis 
is power-spectrum based, the half-ring data sets for the different 
CMB estimators are used. The approach is as follows: 

1. The half-ring temperature maps are multiplied with an ap- 
propriate Galactic and point source mask. 
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2. The cross power spectrum between the t wo halves is esti- 
mated locally using the MASTER approach ( Hivon et al. 2002) 
for 768 discs of diameter 22.5° degrees centered on the pixel 
centers of all the pixels of an N^q = 8 HEALPix map. 

3. We apply the same procedure to the set of 500 simulated 
maps of CMB and noise. 

4. The 768 local spect ra are binned in to blocks of about 100 
multipoles as in Hanse n et al.| ( |2009] ). There are not exactly 
100 multipoles in each block, as the spectra have been ini- 
tially estimated in 16-multipole blocks. 

5. For each 100-multipole block and each disc, the mean power 
from simulations is subtracted and the result is divided by 
the standard deviation. Dividing the spectra by the local stan- 
dard deviation avoids the problem that directions close to the 
Galaxy, where the Galactic mask increases the variance, can 
dominate the statistics due to large fluctuations. 

6. Each 100-multipole block now has an associated map at 
Nside = 8, where each pixel corresponds to the normalised 
power spectrum estimated on the disc centered on that direc- 
tion. 

7. Spherical harmonic transforms are computed for each of the 
maps in order to obtain the dipoles and the dipole direc- 
tions (#,0) of the power asymmetry for each 100-multipole 
block. The alignment of this direction between the differ- 
ent multipole blocks is then a measure of the power spec- 
trum asymmetry. Despite the mask-induced correlations be- 
tween adjacent multipoles, the power spectra estimated in 
100-multipole blocks should be independent and the dipole 
directions of an isotropic field should be random. 

In order to assess the significance of the asymmetry, one has 
to find out whether the distribution of dipole directions for the 
different scales are as random and independent as in the simu- 
lated maps. For this purpose, we define a dispersion angle, # mea n, 
which is the mean angle between all possible combinations of 
100-multipole dipole directions up to a given £ max . We calculate 
#mean(dax = 1500) for the data and compare it to the simula- 
tions. 



Table 22 presents a summary of the power asymmetry results 
from the Planck data processed by the four foreground clean- 
ing methods — Commander -Ruler, NILC, SEVEM and SMICA — 
computed on the U73 mask. To illustrate the effect of the mask in 
the analysis, we also show the result obtained using the SMICA 
data with a smaller mask with / sky = 88% (the CS-SMICA89 
mask which corresponds to the c onfidence mask for that method 
- see Planck Collaboration XII 2013). For comparison, we have 
also included the latest WMAP 9-year result computed with their 
KQ85 mask. 

From Table [22] we see that the result from the mask with 
larger sky coverage, the SMICA data with the CS-SMICA89 
mask, has the highest power asymmetry, the data dispersion an- 
gle of the first 15 100-multipole dipole directions being lower 
than all the 500 simulations. The significance decreases to about 
99.2% confidence level, however, for the U73 mask with a 
smaller sky coverage, except for the case of Commander -Ruler, 
which has an even lower significance. Moreover, the dispersion 
angles among the first 15 100-multipole dipole directions for the 
four methods are consistent. The comparison between the simu- 
lations and the data dispersion angles up to £ max = 1500 is shown 



in Fig. 26 



In Fig.|27]we show the dipole directions of the 15 initial 100- 
multipole bins for the SMICA map with the CS-SMICA89 mask, 
as well as the 6 first 100-multipole bins for WMAP9 data with 
the KQ85 mask (squares). We see that the direction of the first 



Table 22. Summary of the power dipole directions on the sky, up 
to Amx = 1500, as determined from maps of the power spectrum 
estimated from 768 22.5° radius discs and averaged over A£ = 
100 bins. The significance of the power asymmetry, shown in 
the last column, is quantified by the fraction of simulations that 
have smaller clustering of the dipole directions than the data. For 
the Planck analysis we used the 500 FFP6 simulations, while for 
WMAP we used 10000 Gaussian simulations. 



Method 


Mask 


(lb) [°] 


^mean C ] 


1 Idt. (? mean ^. (7 mean 


C-R 


U73 


(231,-5) 


67.8 


11/500 


NILC 


U73 


(223 -1) 


66.1 


4/500 


SEVEM 


U73 


(224,2) 


66.6 


4/500 


SMICA 


U73 


(225,1) 


66.2 


4/500 




WMAP9 


WMAP9 KQ85 


(226,-27) 


33.2 


27/10000 


SMICA 


CS-SMICA89 


(224,0) 


55.8 


0/500 




— C-R 

— NILC 

— SEVEM 
SMICA 

— SMICA** 




100 



200 300 
Simulation Number 



400 



500 



Fig. 26. Dispersion angles of the power spectra dipole direc- 
tions, the mean of the differences of the dipole direction angles, 
up to ^ max = 1500, of the 500 FFP6 simulations compared to 
the Planck data with different foreground cleaning methods. All 
analyses, except SMICA**, are performed with the U73 mask. 
The SMICA** case is for SMICA data with the CS-SMICA89 
mask. 



6 dipoles are similar to the directions found in the WMAP data. 
The preferred directions for WMAP9 and Planck over the range 
£ = 2- 600 are indicated, together with the Planck direction for 
the total range I - 2 - 1500. Finally , the direction of the dipole 
modulation described in Sect. l5.5.2l is also included. Similar be- 
haviour is seen for all of the Planck foreground cleaned maps 
and for the U73 mask, although the scatter between the dipole 
directions increases with increasing sky cut. 

It should be apparent that the asymmetry direction from the 
largest to the smallest angular scales are in general tightly clus- 
tered around the same direction as found for WMAP. However, 
with the Planck data a second preferred direction is also seen 
which is aligned with the CMB dipole direction. This re- 
sult is consis tent with the findings of |Planck Collaboration 
XXVII ( |2013| ), who reports a statistically significant detection 
of Doppler boosting aligned with the CMB dipole at small an- 
gular scales. 
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Fig. 27. Dipole directions for 100-multipole bins of the local 
power spectrum distribution from £ = 2 - 1500 in the SMICA 
map with the CS-SMICA89 mask applied. We also show the to- 
tal direction for £ max = 600 for WMAP9 (black X) and SMICA 
(white X) as well as for £ max = 1500 for SMICA (white big +). 
The stars with different colors correspond to C-R (green), NILC 
(deepskyblue), SEVEM (red) and SMICA (orange) with th e U73 
mask. The best fit dipole modulation direction from Sect. 5.5.2 
is indicated by the white open circle. 



In Fig. 28 we show the Q computed in discs of diameter 
90° centered on the preferred asymmetry dipole direction for £ - 
2-1500 as well as the opposite direction. We can clearly see that 
one spectrum lies systematically above the other over the full 
multipole range, but in particular for the lowest multipoles. Such 
an asymmetry is not seen at the same level of significance when 
the spectra are computed for discs centred on the cosmological 
dipole direction. 

5.5.2. Dipole modulation 

In Sect. |5.5.1| it was shown that the previously reported power 
asymmetry is visible at all multipoles probed by Planck with 
a fairly consistent preferred axis across angular scales. No ex- 
plicit parametric model was assumed in the analysis. In this 
section, however, we revisit the phenomenological model due 
to Gordon et al.| p005| ) considering only large angular scales, 



who proposed that the power asymmetry could be described in 
terms of a multiplicative dipole modulation model of the form 
d = (1 + A p • w)Si S0 + n = Msi S o + n, where A is the dipole ampli- 
tude, p is the dipole direction, n denotes instrumental noise, and 
Si S o is an underlying isotropic CMB field. Both Si S0 and n are as- 
sumed to be Gaussian random fields with covariance matrices S 
and N, respectively. Since Si S0 is assumed to be isotropic, its co- 
variance may be fully specified by some angular power spectrum 

Q,iso- 

In the following we present the results from a direct like- 
lihood analysis of thi s model, similar to t hose described by 
|Eriksen et all J2007a) ; |Hoftuft et all ( [2009] ) for the 3- and 5- 
year WMAP data, respectively. Since this method requires ma- 
trix inversions and determinant evaluations, the computational 
expense scales as 0(N P i x ), and it is therefore feasible only at 
low resolutions. Specifically, we consider maps downgraded to 
a HEALPix pixel resolution of N V [ X = 32, smoothed to angular 
resolutions ranging from 5 to 10°, ensuring sufficient bandwidth 
limitation at this pixelization. All four Planck CMB solutions are 
included in the analysis; however, note that the Galactic plane 
is handled differently in the four approaches. Specifically, for 
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Fig. 28. Top: The power spectra calculated on discs with diam- 
eter 90° for the range £ = 2 - 1500 in the direction of maximal 
asymmetry and its opposite. Bottom: The equivalent plot for the 
direction defined by the comological dipole. The lower panels 
indicate the normaised difference of the spectra from opposing 
directions. 



the Commander map the region inside the corresponding anal- 
ysis mask has been replaced with a Gaussian constrained real- 
ization, eliminating the possibility of bright Galactic residuals to 



leak outside the mask during degradation (Planck Collaboration 



XV 2013 ); for SMICA and NILC a smaller region is replaced with 
a Wiener filter; while for SEVEM no special precautions are taken. 

After degrading each map to the appropriate resolution, we 
add random uniform Gaussian noise of 1//K rms to each pixel to 
regularize the covariance matrix. All pixels inside the U73 mask 
are excluded, and we adopt the difference maps between the raw 
Planck LFI 30 GHz and HFI 353 GHz maps and the SMICA CMB 
solution as two foreground templates, tracing low- and high- 
frequency foregrounds, respecively. We marginalize over these 
Galactic foreground templates, /, as well as four monopole and 
dipole templates, by adding corresonding term of the form aff T 
to the total data covariance matrix, where a is set to a numeri- 
cally large value. 

Before writing down the likelihood for A and /?, a choice 
has to made for the power spectrum, Q s i S0 . We follow [Eriksen| 
et al. ( 2007a| ), and adopt a simple two-parameter amplitude-tilt 

parameter model on the form Q ? i S0 = q (^/^pivot) Q,fid f° r this 
purpose, where the fiducial spectrum, C^m, is the best-fit Planck 
spectrum. The full model therefore includes five free parameters, 
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Table 23. Summary of dipole modulation likelihood results as a function of scale for all four Planck CMB solutions. 



Data set FWHM [°] A (l,b) [°] AlnX Significance 

Commander 5 0.078^ (227,-15) + 19 O 3^ 

NILC 5 0.069^0 (226, -16) + 22 7.1 3.0a 

SEVEM 5 0.066^ J (227, -16) + 24 6.7 2.9cr 

SMICA 5 0.065+J021 (226, -17) + 24 6.6 2.9cr 

WMAP5 ILC 4.5 0.072 + 0.022 (224, -22) + 24 7.3 3.3a 

Commander 6 (223, -16) + 25 6.4 2.8a 

NILC 6 0.062+™*; (223, -19) + 38 4.7 2.3cr 

SEVEM 6 0.060+^5 (225, -19) + 40 4.6 2.2a 

SMICA 6 0.058+°;°^ (223, -21) + 43 4.2 2.1a 

Commander 7 0.062+g^ (223, -8) + 45 4.0 2.0a 

NILC 7 0.055^ (225, -10) + 53 3.4 1.7a 

SEVEM 7 0.055^ (226, -10) + 54 3.3 Ua 

SMICA 7 (226, -11) + 58 2.8 1.5a 

Commander 8 0.(M3+°;°g (218, -15) + 62 2.1 1.2a 

NILC 8 0.049^2 (223, -16) + 59 2.5 1.4a 

SEVEM 8 0.050^2 (223, -15) + 60 2.5 1.4a 

SMICA 8 0.041+™® (225, -16) + 63 2.0 1.1a 

Commander 9 0.068;^ (210, -24) + 52 3.3 1.7a 

NILC 9 0-076 + ™^ (216, -25) + 45 3.9 1.9a 

SEVEM 9 0-078 + ™35 (215, -24) + 43 4.0 2.0a 

SMICA 9 0.070 + ™^ (216, -25) + 50 3.4 1.8a 

WMAP3 ILC 9 0.114 (225,-27) 6.1 2.8a 

Commander 10 0.092^ (215, -29) + 38 4.5 2.2a 

NILC 10 0.098 + ™g (217, -29) + 33 5.0 2.3a 

SEVEM 10 0-103 + ™37 (217, -28) + 30 5.4 2.5a 

SMICA 10 0.094+°;°^ (218, -29) + 37 4.6 2.2a 



namely three dipole parameters and two power spectrum param- 
eters. 

Taking advantage of the fact that both the signal and noise 
are assumed Gaussian, the exact likelhiood may be written down 
in a convenient closed form, 



£(A,p, q, n) oc 



e -\d T (M T SM+N+a^ififJ) l d 

^IMTSM + N + aS//^! 



(39) 



This expression forms the basis of all calculations presented in 
the following. 

Due to the high computational expense associated with these 
evaluations, we do not compute the full joint five-parameter 
model in this analysis, only conditionals of it. However, we 
iterate once in a Gibbs-sampling like approach, by maximiz- 
ing each conditional to obtain an approximation to the full 
maximum- likelihood solution. That is, we first map out the 
dipole likelihood for the 5° FWHM case, fixing the power spec- 
trum at the fiducial spectrum, £(A,p\q = l,n - 0), and locate 
the maximum-likelihood dipole parameters. Then we map out 
the corresponding power spectrum conditional, £(q, n\A m £, p m i). 
Finally, we update the dipole likelihood with these power spec- 
trum parameters, and evaluate the final results. Note that the 



power spectrum and dipole modulation parameters are only 
weakly correlated, and this procedure is therefore close to op- 
timal. Further, the approach is also conservative, in the sense 
that it will always underestimate the significance of the dipole 
modulation model; the derived maximum-likelihood value will 
always lie slightly below the true maximum-likelihood point. 

The results from these calculations are summarized in 
Table [231 listing results for all four Planck CMB maps at an- 



gular scales between 5 and 10° FWHM. For easy reference, we 
also list the results from the corrsponding 3- and 5 -year WMAP 
analyses ( [Eriksen et al.|2007a[|Hoftuft et al.|2009) . Note that the 
former was performed at a HEALPix resolution of Af^e = 16 and 
the latter at an angular resolution of 4.5° FWHM. 



Fig. 30 shows marginals for A, q and n, as derived from the 
Commander CMB solution for all smoothing scales. At least two 
interesting points can be seen here. First, while there is clearly 
significant scatter in the derived dipole modulation amplitude for 
different smoothing scales, as originally pointed out by (Hanson] 
& Lewis | ( |2009| ), all curves appear to be consistent with a single 



value of A ~ 0.07. No other single value fits all scales equally 
well. Second, it is interesting to note that the low-£ power spec- 
trum derived here is consistent, but not without some tension, 
with the fiducial spectrum, (q, n) - (1,0), around 1.5 - 2<r. In 



28 



Planck Collaboration: Isotropy and statistics 




0.05 0.10 0.15 

Modulation amplitude, A 



T3 



1 


1 1 


5 deg 






6 deg 






7 deg 






8 deg - 






9 deg 






10 deg 


^ — 1 


1 1 





0.85 0.90 0.95 1.00 

Power spectrum amplitude, q 




0.0 0.1 

Power spectrum tilt, n 

Fig. 29. Marginal dipole modulation amplitude (top), power 
spectrum amplitude (middle) and power spectrum tilt (bottom) 
probability distributions as a function of smoothing scale, shown 
for the Commander CMB solution. 



particular, there appears to be a slight trend toward a steeper 
and positive spectral index as more weight is put on the larger 
scales, a result already noted by COBE-DMR. The same conclu- 
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Fig. 30. Consistency between component separation algorithms 
as measured by the dipole modulation likelihood. The top 
panel shows the marginal power spectrum amplitude for the 5° 
smoothing scale, the middle panel shows dipole modulation am- 
plitude, and the bottom panel shows the preferred dipole direc- 
tions. The coloured area indicates the 95% confidence region for 
the Commander solution, while the dots shows the maximum- 
posterior directions for the other codes. 



sion is reached using the low-£ Planck likelihood, as described 
in |Planck Collaboration XV| ( [20T3] ). 

In Fig. 30 we compare the results from all four CMB solu- 
tions for the 5° FWHM smoothing scale. Clearly the results are 
consistent despite the use of different algorithms and different 
treatments of the Galactic plane, demonstrating robustness with 
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Fig. 31. Log-likelihood difference between the best-fit dipole 
modulation model and the fiducial isotropic model as a function 
of smoothing scale. Horizontal dashed lines indicate 1, 2 and 3cr 
thresholds. 



respect to the details of the analysis methods. Further, we also 
note that these results are consistent with those derived from the 



5 -year WMAP ILC map by Erikse net al.| ( |2007a] ), demonstrating 
robustness across experiments. On the other hand, it is notable 
that a higher dipole amplitude was found at 9° FWHM for the 3- 
year WMAP ILC map than is observed here, using a larger mask. 

In Fig. [31] we show the log-likelihood difference between 
the derived maximum-likelihood point and the isotropic model, 
A = 0, as a function of smoothing scale. The power spectrum 
parameters are kept fixed at the best-fit values for both points, 
leaving three additional parameters for the dipole model. The 
dashed horizontal lines indicate the 1, 2 and 3<x confidence re- 
gions for three degrees of freedom. As has been noted previously 
in the literature, these significances vary with smoothing scale. 
Taken at face value, the results presented here are suggestive but 
clearly not decisive, resulting in an unchanged situatio with re- 
spect to earlier reports. This is of course not unexpected given 
that WMAP is already strongly cosmic variance limited at these 
angular scales. 

The critical question is whether the trend seen at smaller an- 
gular scales in Fig. 31 continues, or if the apparent likelihood 
peak a t 5° FW HM hap pens to be a local ma ximum. Hanson & 
[Lewis (2009), and later Bennet fet al.| ( |20li) , address this ques- 
tion through a computationally cheaper quadratic estimator, al- 
lowing them to extend a similar analysis to small scales. In doing 
so, they claim that the apparent likelihood peak at 5° is indeed 
a local maximum, and the evidence for the modulation model 
falls off when more data are included. In this respect, it should 
be noted that the dipole modulation model was originally pro- 
posed by Gord on et al.| ( |2005| ) as a simple phenomenological 
characterization of the more general power asymmetry. In par- 
ticular, it assumes that the modulation amplitude, A, is equally 
strong on all scal es. From both the result s shown in Sect. |5.5.1| 
and presented by Hanson & Lewis (2009); Bennett ~et al.| ( |201 1| 7 
this appears not to hold, as the fractional hemispherical power 
difference is clearly smaller at £ > 300 than at £ < 100. On 
the other hand, the preferred directions derived from the current 
\osn-£ analysis is remarkably consistent with the high-£ direction 
derived in Sect. |5.5.l) A proper modulation model may therefore 
need additional spatial structure beyond the simple dipole pro- 



posed by Gordon et al. ( 2005|), as already suggested by Hoftuft 
[eTaLl ( [2009F and |Moss et al.| ( |2011| ). 



5.6. Generalized modulation 

In this section, we study a generalization of the dipolar modu- 
lation field analysed in section [53^2] using the Bipolar Spherical 
Harmonic (BipoSH) formalism. For a statistically isotropy sky, 
the spherical harmonic space two-point correlation matrix is 
diagonal, and, given by the angular power spectrum Q. The 
BipoSH representation provides a natural, mathematically com- 
plete, generalization of the angular power spectrum that captures 
statistical isotropy violations via coefficients that are a com- 
pletely equivalent representation of the spherical harmonic cor- 
relation matrix, 



A LM 



\r LM 



(40) 



This relationship combines the off-diagonal spherical harmonic 
correlations into a bipolar multipole L, M - analogous to the total 
angular momentum addition of states. The CMB angular power 
spectrum corresponds to the L = BipoSH coefficients Q = 

(-i/a^^/V27TT. 

A simple model that results in the violation of statistical 
isotropy arises from the modulation of the of the CMB sky, 



T(n) = T (n)(l + M(n)) , 



(41) 



where T(n) represents the modulated CMB sky, T (n) is the un- 
derlying statistically isotropic random CMB sky and M(n) is a 
fixed, zero-mean, dimensionless, modulation field. The modula- 
tion signal, if any, is expected to be weak and allows quadratic 
terms in M to be neglected. The BipoSH coefficients for the 
modulated CMB field (L > 0) are then given by the following 
expression, 



V 2 



= A^ 2 + m LM G L 



+ Q 2 n^n^ 2 



V47T 



->L0 



(42) 



where A^J? corresponds to the BipoSH coefficients of the un- 
known, but statistically isotropic, unmodulated CMB field, mlm 
are the spherical harmonic coefficients of the modulating field 
(L > 0), Q is the best-fit CMB angular power spectrum and 
= V2£ + 1. The statistically isotropic nature of the unmodu- 
lated CMB sky implies that the expectation values of A^f van- 
ish for (L > 0), leading to the estimator for the modulation field 
harmonics, 



[ { [ 2 



w, 



XLM 

°v 2 



(43) 



denoted by the overhat (Hanson & Lewis 2009). The weights 



wy p for a minimum variance estimate for the modulation field 



correspond to 



Kh = N 



( G L 

(J A LM 



(44) 



Wt t = l.The 



where N is a normalisation chosen such that Y*i x & 
BipoSH representation further allows an estimate of the modu- 
lation field over specific angular scales by windowing regions 
in multipole space in the sum over multipoles £\ , £2 in eqn. HT 



30 



Planck Collaboration: Isotropy and statistics 




Fig. 32. The significance of the modulation power, L(L + 
l)m L /2n, at bipolar multipoles L. The modulation spectra ob- 
tained from the four component separation maps (C-R, NILC, 
SEVEM and SMICA) are consistent with each other. Dipole (L = 1) 
modulation power is detected in all the spectra at a significance 
ranging from 3.7 to 2.9<x. The solid black lines denote the 3<x 
significance thresholds. There is no significant power detected 
at higher multipole of the modulation field 1 < L < 32. 



Fig. 33. The CMB multipole dependence of the BipoSH (mod- 
ulation) power L(L + l)mL/27r can be dissected into bins in £- 
space. This figure plots the measured dipole modulation (L = 1) 
power in CMB multipole bins. We establish that significant 
power in the dipole modulation is limited to £ <e (2, 64) and 
does not extend to the higher CMB multipoles, t, considered. 
The vertical grid lines denote the CMB multipole /'-bins. 



Table 24. This table lists the amplitude and direction of the 
dipole modulation in Galactic coordinates. The measured values 
of the dipole amplitude and direction are consistent for all maps. 
The corresponding dipole power for the SMICA map is seen at a 



detection significance of 3.7<x as shown in Fig. 32 



Map 


Dipole Amplitude 
A 


ihb) [°] 
(cr z = 15.4, a b = 15.1) 


C-R 




(218.9, -21.4) 


NILC 


070 +0 01 


(220.3 , -20.2) 


SEVEM . . 


0065 +o.oii 

U,UUJ -0.011 


(221.7,-21.4) 


SMICA . . 


073 +001 


(217.5 , -20.2) 



This additional information could be very useful in identifying 
the origin of the statistical isotropy violation, which co uld be ei- 
ther cosmological or due to systematic artef acts (see Hajian & 
Souradeep 2003; Hajian & Souradeep 2006). 

First, we limit our analysis to the four low resolution A^e = 
32 CMB maps used in Sect |5.5.2] and reconstruct the modulation 
maps for each of them at the same low resolution. The U73 mask 
is applied to the reconstructed modulation maps before comput- 
ing m LM - The pseudo-power m L is corrected for the mask applied 
to the modulation maps. Specifically for the case of dipole mod- 
ulation, the pseudo -power ml is related to the dipole amplitude 
by A = 1.5 ^mi/n. 

A dipole modulation (L = 1) signal is detected at 3<x sig- 
nificance in all the maps, as shown in Fig. 32 The amplitude 
and direction of the dipole m odulati on match those obtained via 
a likelihood analysis in Sect. |5.5.2| The BipoSH representation 



of modulation confirms the dipole modulation signal found in 
the low-resolution map. Since this approach allows the recon- 
struction of any general small amplitude modulation field, the 
BipoSH representation places constraints on the power in the 
modulation field at all higher (bipolar) multipoles allowed by 
the resolution of the CMB maps. 

We then extend the analysis to higher resolution using maps 
at A^ S i de = 256 for Commander and Af side = 2048 for NILC, SEVEM 
and SMICA in order to study the above effect in more detail. 
We repeat the analysis on these higher resolution maps using 
the U73 mask. Contrary to our expectations based on a scale- 
independent (i.e., no /'-dependence) model, the significance of 
the dipole does not increase in the high resolution maps. We then 
subdivide the Grange up to £ max = 384 into uniform bins of size 
A£ = 64. As seen in Fig. 33 we recover the dipole modulation at 
over ~ 3<t significance only for the lowest bin I e (2, 64). This 
is consistent with the results in Sect. |5. 5. 2] and the BipoSH anal- 
ysis on the corresponding low resolution maps shown in Fig. [32] 
However, the amplitude of the dipole is consistent with zero 
within 3<r for all of the higher ^-bins considered. This suggests 
that the simple modulation model in Eqn. 41 is inadequate and 
should minimally allow for the amplitude, A(^), of the dipole to 
depend on CMB multipole, I. Although this may appear to be 
a more complex model, it does not necessarily lack motivation. 
It is readily conceivable that physical mechanisms that cause a 
dipolar modulation of the random CMB sky would be scale de- 
pendent and possibly significant only at low wavenumbers. More 
importantly, such a dip ole modulation has also been noted in low 
resolution WM AP data flEriksen et al.|2007a[|Hoftuft et al.|2009] ). 
More recently, Bennett et al. ( |2011| ) also comment (without be- 
ing quantitative) that the effect is present in the WMAP maps 
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A LM 



but limited to low £ and conclude that the £ dependence rules 
out a simple modulation explanation. The fact that two indepen- 
dent experiments find this intriguing statistical isotropy violation 
points to a non-instrumental origin. 

It is, of course, possible to extract the BipoSH coefficients 
, up to the maximum multipole £ m2LX allowed by the full res- 
olution Planck maps at modest computational expense. This al- 
lows us to address a specific indication of s tatistical isotropy vi- 
olation previously reported in the literature. Bennett et al. ( |2011| ) 
found nonzero BipoSH power spectra, A® and A^ +2 at very 
high statistical significance in the WMAP maps as determined 
in ecliptic coordinates, corresponding to a quadrupolar power 
asymmetry in the CMB sky. The BipoSH spectra peaked at £ ~ 
250, and the differences in the BipoSH signal determined from 
two different frequency bands indicated a non-cosmological ori- 
gin. Furthermore, the azimuthal symmetry of this BipoSH signal 
in ecliptic coordinates suggested that it had its o rigin in some 
unaccounted-for systemat ic effect. The findings of Han son et al. 
2010| ); Joshi et aL| ( 2012| ) strongly suggest that the signal arises 



due to an incomplete treatment of beam asymmetries in the data. 
Bennett et al. ( 2012 ) have subsequently noted that analysis of the 
WMAP9 beam-deconvolved maps no longer detects the signal. 

We have computed the A 2 ® and A^ +2 in Ecliptic coordinates 
for the full resolution Planck CMB maps as shown in Fig. 34 
The analysis yields no evidence for BipoSH coefficients that de- 
viate significantly from zero. This provides conclusive observa- 
tional evidence from independent CMB measurements that the 
WMAP result could have only arisen due to instrumental arte- 
facts in that data set. 



5.7. Parity asymmetry 

5.7.1 . Point-parity asymmetry 

The CMB sky map may be considered as the sum of even and 
odd parity functions. Previously, an odd point-parity preference 
(hereafter parity asymmetry) was observed in the WMAP 3-, 
5- and 7-year data releases ([Land & Magueijo 2005 b[ [Kim & 
| Naselsky|2010a]|Naselsky et al.|2012[|Kim & Naselsky|2010b| 
Gruppus o et al.||2011| ). In this section we investigate the parity 
asymmetry for the Planck temperature anisotropy power spectra 
derived with a quadratic maximum likelihood (QML) estima- 
tor applied to the Commander -Ruler, NILC, SEVEM, and SMICA 
maps at N S ^ Q = 32, and with a pseudo-Q estimator at N S ^ Q = 64. 

From the CMB anisotropy field defined on the sky, T(n), we 
may construct symmetric and antisymmetric functions using the 
coordinate inversion n — > -n: 



T + (n) : 



T(n) + T(-n) 



T~(n) 



T(n) - T(-n) 



(45) 



Therefore, T + (n) and T~(n) have even and odd parity, respec- 
tively. When combined with the parity property of spherical har- 
monics, Yi m (ri) - (-1/ Y{ m (-ri), we obtain: 



T + (n) = Yj a ^Y, m (n)Y\£\ 

l,m 

T~{n) = ^a im Y €m {n)Y~(£), 



(46) 



e,m 



where n is an integer, and Y + (£) = cos 2 (y), and Y (£) = 
sin 2 (f). 

A significant power asymmetry between even and odd mul- 
tipoles may thus be interpreted as a preference for a particular 
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Fig. 34. The BipoSH power spectra Af t and Af^ +2 obtained 
from the four component separation maps (C-R, NILC, SEVEM 
& SMICA) are consistent with each other. Note that no signifi- 
cant (> 3<x) detections are found. This independently establishes 
the fact that the quadrupolar BipoSH detections made by WMAP 
were due to specific instrument systematics. 



parity of the anisotropy pattern, connected to the parity asymme- 
try of the metric perturbations at scales above 1-4 Gpc ( |Kim & 
|Naselsky|2010a| ). For investigation of the parity asymmetry w e 
may consider the following quantities ( Kim & Naselsky 2010a): 



p-(0 



8d) = 



n=2 

Z r ~ (n) 

n=2 

P~(£) 



In 



n{n + 1) 
In 



Cn-> 



(47) 



where P + and P~ are the sum of n(n + l)/2n C n for even and 
odd multipoles respectively; the ratio P + /P~ is associated with 
the degree of parity asymmetry, where a value of P + /P~ < 1 
indicates odd-parity preference, and P + /P~ > 1 indicates even- 
parity preference. 

Following ( |Kim & N aselsky 2010a), we will discuss the 
range of multipoles 2 < £ < 30, which belongs to the Sachs- 
Wolfe plateu of the TT power spectrum, where £(£ + 1)Q ~ 
const. In order to make a rigorous assessment of the statistical 
significance of parity asymmetry at low £, we have compared 
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Fig. 35. Top panel: the parity estimator g(fi) versus t for 
Commander -Ruler (black diamonds), NILC (red diamonds), 
SEVEM (blue diamonds) and SMICA (green diamonds). Bottom 
panel: /?-value for C-R (black solid line), NILC (red line), SEVEM 
(blue line), and SMICA (green line). 



g{€) for the Planck power spectra with 10 3 simulated CMB maps 
based on the fiducial Planck cosmolo gical model. We comp ute 
power spectra using a QML-estimator ( Gruppu so et al.|2009] ) as 
applied to data at A^e = 32 with the U73 mask applied. This 
yields practically identical power spectrum results for the same 
Grange determined with a pseudo-Q estimator applied to maps 
at Af S ide = 64. 



In Fig 35 we show the g(£)-parameter for the Planck power 
spectra and the corresponding p- values. The p- value denotes the 
fraction of simulations in which the obtained value of P + /P~ 
is as low as that observed in the data. Note that the results 
from the different Planck CMB maps yield consistent shapes for 
the g(£) and /^-parameters. The parity asymmetry at I - 22 
is most anomalous, with a corresponding /7-value in the range 
0.002 — 0.004). Finally, the statistical significance of the parity 
asymmetry (i.e., low p- value) increases when we increase ^ max 
up to 22-25. Therefore, the odd parity preference cannot simply 
be attributed to the low quadrupole power. It is plausible the low 
quadrupole power is not an isolated anomaly, but that it shares 
an origin with the odd parity preference (see for details (|Kim &| 
Naselsky|201Qa]|Naselsky et al.|2012HKim & Naselsky|2010b| )). 



5.7.2. Mirror Parity 

In this section we investigate the properties of the Planck tem- 
perature low-resolution maps under reflection with respect to 
a plane. This search for hidden mirror symmetries and anti- 
symmetries co mplem ent s the tests for pa rity asymmetry, pre- 
sented in Sect. 5.7.1 Starobinsky ( |1993| ) showed how a hid- 
den mirror symmetry might be connected to the non-compact T l 
topology, or to a compact T 3 topology in which one topological 
scale is much less than the others. The CMB pattern would then 
exhibit a mirror symmetry with respect to the plane defined by 
the two large dimensions. Mirror symmetry has been searched 
for in the COBE-DMR data in |de Oliveira-Costa et all ( (1996] ), 
resulting in a lower limit for the scale of the compact dimen- 
sionas 4 Gpc (see also Gurzad yan et al.|2007[|B en-David et aT] 
2012 for other more recent analysis). Finelli et al. (2012|) anal- 
ysed hidden mirror symmetry and anti-symmetry properties of 
the WMAP 7 -year ILC temperature map, finding a preferred di- 
rection that could be considerede anomalous at the 93 % con- 
fidence level with anti-symmetry properties. This direction lies 
close to the one defining the hemispherical asymmetry. 

Following Finelli et al. ( 2012| ), we consider the following 
estimators: 



1 x^\l/6T ST \ 



vpix . j 



(48) 



where the sum is meant over the observed pixels, Af p i x , ST/T(rij) 
is the CMB temperature anisotropy measured at the pixel pointed 
by the unit vector rij, and n k is the opposite direction of nj with 
respect to the plane defined by w„ i.e. 



n k 



nj - 2 (m nj)fii . 



(49) 



We compute the quantities S ± for each of the 3072 directions 
defined by HEALPix resolution A^de = 16 map, by allowing the 
j and k indices to run over the unmasked pixels of the low reso- 
lution foreground cleaned maps. We perform the same analysis 
on 1000 simulated skies and store the minimum and maximum 
value for each of these. 

The minimum value for the S + estimator is reached for the 
plane defined by Galactic coordinates (#, cp) = (104°, 262°), with 
a significance of 0.8% (Commander -Ruler), 0.5% (NILC), 9.6% 
(SEVEM), and 1.2% (SMICA). The top panel of Fig. [36] shows the 
minimum value of S + for each of the four methods and com- 
pared to the MC simulations computed for Commander -Ruler, 
which is considered to be representative. 

The minimum value for the S ~ estimator is found for a direc- 
tion close to that associated with the cosmological dipole. It is 
not statistically significant for any of the CMB maps (see bottom 
panel of Fig. ([36])). 

The anomalous anti- symmetry direction found in the Planck 
CMB data is close to that found for the dipolar modulation in 
Sect. |5.5| suggesting some connection between them. The di- 
rections which minimize S + and S~ for Planck are the same 



as those found for the WMAP 7-year ILC map in |Finelli et aL 
( |20T2] ). 



5.8. The Cold Spot 



The Cold Spo t was identified in the WMAP first year data ( Vielva 
|et al. 2004) through the estimation of the kurtosi s of the 



Martinez- 



Spher ical Mexican H at Wavelet (SMHW, e.g. 
Gonzalez et al. 2002) coefficients, and confirmed ( [Cruz et aL 
2005 1) by analysing the area of the SMHW coefficients 
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Fig. 36. Top panel: the S + statistic. The vertical lines show the 
minimum value for the estimator as computed on low resolu- 
tion C-R, NILC, SEVEM, and SMICA maps. The grey histogram 
shows the same quantity computed from 1000 simulated maps 
processed by C-R. Bottom panel: as above for S~. 



Table 25. Upper tail probability (UTP, in %) associated to the 
cold (left) and hot (right) areas. Results are given for a v > 
4<tr threshold and for the four Planck CMB maps. The three 
most significant scales associated to the Cold Spot are shown. 
Analysis performed on the exclusions masks associated with the 
U73 mask. 





Scale (R) 


C-R 


NILC 


SEVEM 


SMICA 




['] 


UTP (%) UTP (%) UTP (%) UTP (%) 


cold area 


200 


1.6 


1.1 


1.2 


1.1 




250 


0.3 


0.3 


0.3 


0.3 




300 


0.3 


0.3 


0.3 


0.3 


hot area 


200 


2.3 


1.6 


1.8 


1.6 




250 


2.7 


2.2 


2.4 


2.2 




300 


4.9 


3.7 


4.1 


3.8 



Table 26. Upper tail probability (in %) associated to the cold 
(left) and hot (right) areas. Results are given for a v > 4cr R 
threshold and for the four Planck CMB maps. The three 
most significant scales associated to the Cold Spot are shown. 
Analysis performed on the exclusions masks associated with the 
G70 mask. N/A indicates that no area above that threshold was 
found on the data. 





Scale (R) 


C-R 


NILC 


SEVEM 


SMICA 




['] 


UTP (%) UTP (%) UTP (%) UTP (%) 


cold area 


200 


1.0 


1.0 


0.9 


1.0 




250 


0.3 


0.3 


0.3 


0.3 




300 


0.2 


0.2 


0.3 


0.2 


hot area 


200 


15.1 


14.5 


14.6 


14.5 




250 


N/A 


N/A 


N/A 


N/A 




300 


N/A 


N/A 


N/A 


N/A 



above/below a given threshold. Since its detection, the Cold Spot 
has been extensively st udied and verified with a large battery of 
statistical probes (e.g., Mukherjee & Wang 2004; |Cay6n et al.| 
2005; McEw en et al.|2005HCruz et al.|2007al|Rath et al.|2007b 



Vielva et al.|2007t |Pietrobon et al.|2008[ |Gurzadyan et al.|2009 
Rossmanith et al. 200 9b|). A completereview of the Cold Spot 
can be found in|Vielva (2010), including a discussion on possible 
explanations of its nature. 

The analysis of the kurtosis of the SMHW coefficients has al- 
ready been addressed in Sect. 4.5 We have checked that the kur- 
tosis of the coefficients corresponding to the four Planck cleaned 
frequency maps is larger than the expected value obtained from 
simulations, with a modified upper tail probability of around 
0.01. This is compatible with the value obtained from WMAP. 

Nevertheless, the Cold Spot is more robustly described in 
terms of a morphological quantity: the area of the SMHW co- 
efficients above/below a given threshold. At a given scale R and 
threshold v, the cold (A~ v ) and hot (A+ y ) areas of the SMHW 
coefficients are defined as: 



A~ R y ee #{aj T (R,p)<-v} 
A + R y ee #{co T (R,p)<+v} 



(50) 



where # represents the number operator, i.e, it indicates for how 
many pixels p, the specific condition defined between the braces 
is satisfied. 



Table [25] summarises the results for the hot and cold areas 
determined for the four CMB maps analysed with the U73 mask 
(and its associated exclusions masks). The cold area is anoma- 
lous at scales between R = 200 and R = 300', similar to the 
sizes already highlighted with the kurtosis analysis. We see that 
the higher the threshold, the smaller the upper tail probability as- 
sociated with the Planck CMB map. In particular, the cold area 
has a upper tail probability of 0.003 at v > 4cr R and for R = 300'. 

Notice that the most significant deviation comes from the 
cold area, although the hot area is marginally compatible. 
However, the cold area represents the most robust detection of an 
anomaly, since it is robust to the mask employed (see Tables 26 
and|27). 

The information provided in the previous Tables is also rep- 
resented (for the R = 300' scale) in Fig. 37 In these nine panels 
we show the anomalous cold (in blue) and hot (in red) areas for 
thresholds v > 3.0cr#, v > 3.5cr# and v > 4.0cr# as determined 
from the SMICA map. For the two largest thresholds, the cold 
area corresponds to the Cold Spot, whereas the red area at 3. Oct 
has a lready been identified in the WMAP data (e.g., Vielva et al. 
2007) as an anomalous hot spot. From these analyses it is clear 
that the Cold Spot anomaly is present in both the WMAP and 
Planck data. 

5.9. Interpretation of anomalies 

The results presented here in Sect, [^demonstrate that many fea- 
tures previously observed in the WMAP data are present also in 
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Fig. 37. SMHW coefficients at R = 300 arc minutes, and thresholds of 3.0cr (left), 3.5<x (middle), and 4. Oct (right). Results for the 
three masks considered in the analysis are shown: U73 mask (top), CG70 (middle) and CG60 (bottom). 



Table 27. Upper tail probabilities (in %) associated with the 
cold (left) and hot (right) areas. Results are given for a v > 
4cr R threshold and for the four Planck CMB maps. The three 
most significant scales associated to the Cold Spot are shown. 
Analysis performed on the exclusions masks associated with the 
CG60 mask. N/A indicates that no area above that threshold was 
found on the data. 



Scale (R) 


C-R 


NILC 


SEVEM 


SMICA 


['] 


UTP (%) UTP (%) UTP (%) UTP (%) 


cold area 200 


1.1 


0.9 


0.8 


0.9 


250 


0.1 


0.1 


0.1 


0.1 


300 


0.1 


0.1 


0.1 


0.1 


hot area 200 


N/A 


N/A 


N/A 


N/A 


250 


N/A 


N/A 


N/A 


N/A 


300 


N/A 


N/A 


N/A 


N/A 



the Planck sky. This agreement between two independent exper- 
iments effectively rules out the possibility that their origin lies in 
systematic artefacts present in either data set. In particular, there 
is evidence for a violation of statistical isotropy at least on large 
angular scales in the context of the Planck fiducial sky model. 
Moreover, a dipolar power asymmetry may extend to scales cor- 
responding to £ ^ 1500, whilst fits to a model containing a dipole 
modulation field yield results in excess of 3<x significance. In ad- 
dition, there is evidence from such fits that the low-/' spectrum of 
the Planck data departs from the fiducial spectrum in both am- 
plitude and slope. These results could have profound implica- 



tions for cosmology. It is therefore pertinent to consider whether 
a model can be proposed to provide a common origin for the 
anomalies. 

The microwave sky is manifestly non-Gaussian and 
anisotropic, with known contributions from Galactic astrophys- 
ical foregrounds, lensing of CMB anistropies by the intervening 
matter distribution, and the ISW. However, the excellent per- 
formances of the component separation algorithms used here 
in rejecting diffuse foregrounds argues strongly against known 
Galactic emission as the source of the anomalies. 



Schwarz et aTlj2004) ,[Copi et al.|p007 



), Paris etal. (2011 

and Hansei Tet al.| ( |2012| ) suggested that diffuse Solar System 



emission could contribute to the observed structure on large 
angular scales, although it is not expected that the classical 
Zodiacal Light Emission or Kuiper Belt objects are responsible. 
Planck Colla boration XIV| ( |2013| ) presents the current Planck 
contribution to the modelling of the Zodiacal cloud. 

Another possibility is that the ano malies have their origin 
in the local Universe. According to Francis & Peacock ( 2009 ), 
the removal of the ISW signal originating within the volume at 
z < 0.3 from WMAP data reduces the significance of the appar- 
ent alignment betwe en the CMB quad rupole and octopole and 
the Cold Spot. [Efstathio u et aT] ([2010 ) have used the same cor- 
rection to yield an increase in the structure of the two-point cor- 
relation function for angular separations less than 60°, that had 
been noted as apparently anomalous since the first WMAP data 
release. A future possibility is that Planck itself will be able to 
reconstruct the ISW signal and test its impact on issues related to 
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0.0 



6.0 



Fig. 38. Same as Fig. 24 but with the best fit Bianchi template 



subtracted from the SMICA map. 



isotropy and non-Gaussianity. [Planck Collaboration XIX ( 2013| ) 
presents maps of the effect based on the current data release. 

Of more interest to us is that the anomalies are genuinely 
cosmological in origin. In that context, obvious candidate mod- 
els include those with simply or multi-connected topology. In a 
companion paper (Planck Collaborati on XXVl||2013| ), a subset 
of such models are considered and the signatures of their spe- 
cific correlation structures on the sky are searched for. However, 
no detections are found, but rather the scale of topology is lim- 
ited to be of order the diameter of the last- scattering surface or 
greater. More interestingly, they reconsider Bianchi Vllh mod- 
els that were previously demo nstrated to show statistical cor- 
relation with the WMAP data ([Jaffe et al.|[2QQ5] [2006} |Bridges| 
et al.|2007 McEwen et al.|2013| ), albeit with parameters incon- 
sistent with standard cosmological parameters. In this new anal- 
ysis, the Bianchi parameters are physically coupled to the cos- 
mological ones, yielding no evidence for a Bianchi VII/* cosmol- 
ogy. However, as before, when treated simply as a template for 
structure contained in the CMB sky, a best-fit pattern is found 
to b e in good agreement with the old results. Previous analy- 
ses ( [Jaffe et al.||2005l |Cay6n et al.|2006| |McEwen et aLpOOS] ) 
have shown that when the CMB sky is corrected for such a tem- 
plate, many of the large-scale anomalies are no longer present 
at a statistically significant level. It is likely that such an ef- 
fect will persist for Bianchi-corrected Planck data, and we have 
made an explicit test as to whether the best-fit Bianchi template 
can also explain the presence of phase corr elatio ns. We therefore 
repeated the surrogate analysis from Sect. |5.4| for the appropri- 
ately corrected SMICA map. Figure [38] presents the result for the 
corresponding significance map. It is clear that the signature for 
hemispherical asymmetry is drastically reduced, thereby render- 
ing the signal formally statistically insignificant. Thus, the best- 
fit Bianchi model can also account for the asymmetries induced 
by higher order phase correlations. It should also be noted that 
subtracting the best-fit Bianchi template from the data, outside 
the U73 mask, explains the anomalous skewness and kurtosis 
values but not the variance, for which the corresponding lower 
tail probabilities are 0.008, 0.166, and 0.306, respectively. Given 
the lack of consistency of the physical parameters of the model 
with the Planck cosmological model, the results obtained using 
Bianchi-subtracted input maps might be considered moot, how- 
ever, the morphology of the maps may provide insight into the 
type of underlying structures associated with the anomalies. 

Although the Cold Spot is also rendered statistically insignif- 
icant by the Bianchi template, other possible explanations about 
its nature have been proposed, including the late evolution of 



the large-scale structure (e.g., Inoue & Silk 2006, 2007), the 
Sunyaev-Zeldovich effect (e.g., |Cruz et al.|2008| ), residual fore- 
grounds ( Cruz et al. 2006 ), gravitational lensing ( Das & Spergel 
2009 ), or a cosmic texture (e.g., |Cruz~ et al. 2007b|). 

The presence of primordial magnetic fields (PMFs) due to 
either pre- or post-recombination mechanisms could also pro- 
vide a physical basis for some of the anomalies discussed in 
this paper. Specifically, PMFs with coherence scales compara- 
ble to the present day horizon could result in Alfven waves in 
the early Universe that generate specific signatures on the sky 
via the Doppler and integrated Sachs- Wolfe effects. In particu- 
lar, a preferred angular direction in the CMB anisotropy can be 
induced ( [Durrer et aL][T998l |Kim & NaseTsk^|2009l ) that leads 
to structure in the spheric al harmonic mode correlation matrix 
( Kahniashvili et al. 2008| ). Appendix [A] presents a search for 



the predicted correlations between harmonic modes separated by 
A£ = 0, ±2, and Am - 0, ±2, allowing constraints to be placed 
on the Alfven wave amplitude. Further constraints on PMFs 
ba sed on the power spectrum and bisp ectr um have been provided 
in |Planck Collaboration XVll ( [20T3] ) and |Planck Collaboration! 
|XXIV|p013| ), respectively. 

To conclude, when analysing a data set as complex and rich 
as that provided by Planck, some statistical outliers will be ex- 
pected. However, it should be clear that the evidence for some of 
the large-angular scale anomalies is significant indeed, yet few 
physically compelling models have been proposed to account for 
them, and none so far that provide a common origin. The dipole 
modulation model that was analysed here was phenomenologi- 
cally motivated, and is detected in the data at relatively high sig- 
nificance. Whether it can resolve the anomalous nature of other 
observed features remains to be evaluated. 



6. Implications for Q and cosmological parameter 
estimation 

The approach to Q estimation, the construction of the Planck 
likelihood and subsequent inference of cos mological parameters 
are describ e d in the accompanying pape rs Planck C ollaboration 
[XVl ( |20T3] ); |Planck Collaboration XVllp013| ). For these stud- 
ies, specific assumptions are made about the isotropy and 
Gaussianity of the primordial fluctuations observed in the CMB. 
The latter in particular seems to be well- supported by the com- 
prehensive set of tests applied to the Planck data in Sect. [4] The 
most significant discr epancies are seen in association with the 
Cold Spot (Sect. 5.8), which constitutes a localized region of 
sky of about 10° in size. Its impact on cosmological parameters 
is then likely to be relatively insignificant, and masking of the 
region could easily test this assertion. 

It is well-known that the quadrupole and octopole have low- 
amplitudes relative to the best-fit cosmological power- spectrum. 
The contribution of those multipoles to cosmological parameter 
estimation is very small due to the associated cosmic variance on 
these scales, an d thu s the direct impact of their alignment (as dis- 
cuss ed in Sect. |5.1|) is also likely to be small. Remarkably, how- 
ever, Planck Collaboration XV (2013 ) presents evidence that the 
low-£ multipole range from 2-30 is coherently low, and is not 
well accounted for by the standard ACDM model. Moreover, 
this conclusion is a consequence of the fact that the cosmolog- 
ical parameters are strongly influenced by the i =1000-1500 
range, previously inaccessible to WMAP. Consistent findings 
have been presented here in the form of the low- variance of the 
data in Sects. pTT1 and |5.2| although this is largely driven by the 
quadupole and octopole alone. 
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The question therefore remains as to whether there is a 
deeper connection with the cosmological anomalies seen in both 
the WMAP and Planck data sets particularly on large angu- 
lar scales. Indeed, the hemispherical asy mm etry and dipolar 
power modulation discussed in Sects. 5.3 and |5.5| respectively 
could have a more important impact in that they directly address 
whether a broader class of cosmological models should be con- 
sidered. Indeed, the low-^ signature seen in the data has pre- 
viously been associated with missing power in a Universe with 
simply- or multiply-connected topology. However, there are spe- 
cific morphological signatures of such topologies that have not 
been detected in the Planck data (Planck Collaborati on XXVl| 
|20T3] ). 

However, the phenomenologically motivated dipole modula- 
tion model due to |Gordon et al.| ( [2005]) y ields a significant fit to 
the data, as seen in Sects . [53^21 and [5^6] The former also shows 
some evidence for a departure from the Planck fiducial power 
spectrum in both amplitude and slope. Both of these analyses are 
in good agreement in terms of the direction of the dipolar mod- 
ulation field with the model independent dipolar power modula- 
tion analysis of Sect. |5.5.1| 

A qualitative exploration as to how these features are re- 
flected in the low-^ power spectrum is provided in Fig. 39 
Specifically, the plot presents the angular power spectra com- 
puted using a quadratic maximum-likelihood (QML) estimator 
( [Pad et al.||2QT0l |20T3] ) from the Af side = 16 SMICA map after 
application of the U73 mask used in this paper. The Planck fidu- 
cial power spectrum is also shown for comparison. Clearly, there 
is a deficit of power as expected when no further partitioning of 
the sky is applied. However, further interesting properties of the 
data are revealed when spectra are computed for the two oppos- 



ing hemispheres defined by the preferred direction in Sect. 5.5.1 
In the positive direction, there is improved agreement between 
the derived spectrum and the Planck fiducial sky, but with an in- 
teresting oscillation between odd and even modes. For the neg- 
ative direction, an overall suppression of power is again seen. It 
would be interesting to test the connection between these spec- 
tral features and the phase correlations detected i n Sec t. |5.4| or 
the evidence for parity violation presented in Sect. 5J_ The ob- 
servations may, in part, reflect the presence of visually striking 
features noted by Bennett et al. ( 201 1) — the four elongated cold 
fingers stretching from near the Galactic equator to the south 
Galactic pole and a prominent cold spot near the center of the 
map. 



However, Fig. [28] and the corresponding analysis suggest 
that the asymmetry in power between hemispheres extends to 
much smaller angular scales. Whether such a property of the data 
would have implications for parameter estimation may yet need 
further exploration. 

7. Conclusions 

In this paper, we have tested the statistical isotropy and 
Gaussianity of the CMB using data from the Planck satel- 
lite. We have demonstrated that little evidence is seen for non- 
Gaussianity, although some deviations from isotropy are found. 

Most of the tests performed in Sect. [4] showed an overall 
consistency with the null hypothesis, as represented by a set of 
realistic simulations based on a Planck fiducial sky model and 
including the secondary ISW-lensin g effect (as detected for the 
first time with the Planck data, see |Planck C ollaboration XIX 
2013). However, two important exceptions were seen. The vari- 
ance of the CMB signal was found to be significantly lower than 
expected, with the anomalously low signal seemingly localized 
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Fig. 39. Angular power spectrum at large angular scales com- 
puted on opposing hemispheres defined by the maximal asym- 
metry. 



in the northern ecliptic hemisphere (Sect. |4.1| ). This result was 
also confirmed with the low variance of the wavelet coefficients 



that was seen on scales above a few degrees (see Sect. 4.5). 
Moreover, a significant deviation from Gaussianity was found 
in the form of a positive kurtosis of the wavelet coefficients. 

These results correspond to statistical features on large angu- 
lar scales where numerous anomalies were previously observed 
in the WMAP data. In Sect. [5j we revisited these in the light of 
the Planck data and found results in excellent agreement with 
those for WMAP. In particular, the most significant anomalies, 
namely the quadrupole-octopole alignment (Sect. |5.1| ), hemi- 
spherical asymmetry (Sect. |5.3| ), pha se correlations (Sect. |5.4| ), 
dipolar power modulation (Sect. |5.5| ), generalized power modu- 
lation (Sect. 5.6 ), parity asymmetry (Sect. 5.7 ) and the Cold Spot 
(Sect. |5.8| ) have been confirmed with the Planck data. Attempts 
to explain the observed features in terms of systematic artefacts, 
local astrophysical sources of emission, or structure in the local 
Universe have not been successful. It is clear that these anoma- 
lies represent real features of the CMB sky. 

However, it is difficult to make a detailed interpretation of the 
anomalies in the absence of theoretical models, in particular with 
regard to the role of a posterior choices. Nevertheless, Planck 
does offer new possibilities to check the a posteriori claims in 
this context as a consequence of its superior multipole content 
that cannot easily been probed by WMAP. This is partic ularly 
relevant for the power asymmetry studies — Sect. 5.5.1 found 
that the same direction was preferred at I > 600 as for I < 50, 
which should mitigate in part such criticisms. 

Phenomenological models have been suggested to account 
for the observations. The d ipolar power modulation app roach 
due to Gordon et al. ( 2005| ) was explicitly tested in Sect. 5.5.2 
and found to represent a good fit to the large scale asymmetry, 
corresponding to a detection at about 3<x significance. This re- 
sult was con firmed by the more generalized modulation study 
in Sect. |5.6| which also ruled out the presence of modulation 
fields of higher order. Alternatively, a Bianchi template f it to 
the data performed in Planck Collaboration XXVI (2013 ) can 
provide a good fit to the hemispherical asymmetry, the Cold 
Spot and the phase correlations, but corresponds to values of 
the cosmological parameters incompatible with those derived in 
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Planck Collaborat ion XVI ( |2013 ). Clearly, these do not provide 
complete and satisfactory explanations for the observations, and 
more physically motivated models should be sought. 

This may also be indicated by the cosmological parameter 



studies presented in Planck Collaboration XV (2013); Planck 
[Coll aboration XVI ( 2013). Here, it was demonstrated that while 
the power-spectrum determined from the Planck temperature 
data is extremely consistent with a basic 6-parameter ACDM 
model, the low-£ multipoles {I < 30) deviate from the best-fit 
model although at a significance that does not appear to exceed 
2.1cr. However, this is precisely the regime where many of the 
anomalies presented in this paper seem to manifest themselves, 
and where qualitatively interesting differences are observed in 
the power- spectra for two hemispheres defined by the preferred 
direction for the dipolar power modulation. 

Finally, it is expected that the polarization data that will be- 
come available with the 2014 data release should provide valu- 
able information on the nature of the CMB anomalies. Then, 
the presence, or even absence, of a specific signature in the data 
should help to elucidate the physical mechanism that is causing 
the anomaly (see|Vielva et al.||2011l Frommert & EnBlin 2010| 
and|Dvorkin et al. 2008 for examples related to the Cold Spot, 



mode alignment, and dipolar modulation, respectively) In par- 
ticular, a deviation of isotropy present at recombination should 
be reflected in both the temperature and polarization data with 
a correlated signal. It may be that the statistical anomalies cur- 
rently described in this paper are a hint of more profound physi- 
cal phenomena that are yet to be revealed. 
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Appendix A: Constraints on Alfven waves 

Observations of synchrotron emission and Faraday rotation pro- 
vide increasing evidence that large-scale astrophysical systems 
in the Universe are pervaded by magnetic fields. These huge sys- 
tems include hy-a forests and intercluster regions (see Kronberg 
2009, for a review). Both pre- and post-recombination mech- 
anisms could result in a background of nano-gauss fields that 
might be detectable in large-scale structures or the CMB, al- 
though at present no imprints of these Primordial Magnetic 
Fields (PMFs) have been detected therein. 

Here, we report our findings based on an analysis of the 
Planck data to search for the the predicted signature of statisti- 
cal anisotropy due to PMFs. Specifically, PMFs with coherence 
scales comparable to the present day horizon may induce and 
sustain Alfven waves in the early Universe that can leave ob- 
servable imprints on the CMB via the Doppler and integrated 
Sachs- Wolfe effects. In particular, this results in a preferred an- 
gular direction in the CMB anisotropy, therefore breaking statis- 
tical isotropy (Durrer et al.|1998||Kim & Naselsky |2009|). 



Durrer et al.|(|1998|) showed that cosmological Alfven waves 



generate a fractional CMB anisotropy for a Fourier mode k: 



ST 



(h, k) « n • n(&, 7/i ast ) = n • £! v A kni ast B • k 



(A.l) 



where n denotes sky direction, B is a unit vector in the direc- 
tion of the coherent PMF, Q(k, 7/i ast ) is the Gauge invariant lin- 
ear combination associated with vector perturbations, rj\ ast de- 
notes the conformal time at the moment of baryon-photon de- 
coupling, and Tq is the CMB monopole temperature 2.7255 K 
(Fixsen 2009]). [Durrer et al.1 fl998) assumed that the vector per- 
turbations are initially created by some random stochastic PMF 
and have the following statistical properties over an ensemble of 
universes: 

(Q! (k) fljj(k)> = (Sij - h kj)P(k). (A.2) 

Here, P(k) is the power spectrum assumed to follow a simple 
power law: 



P(k) = A v 



(A.3) 



where ko is a pivot wavenumber set to 0.05/Mpc in this an alysis. 
The Alfven wave velocity is given by ( [Durrer et al.|1998] ): 

B B 

va = — - ~ 4 x 10" 4 ^_ Q ^ - , (A.4) 



2 ^n(p r + p r ) 



10" 9 Gauss' 
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where p r and p r are the density and the pressure of photons, and 
the speed of light c is set to 1 . 

Kahniash vili et al. ( 2008| ) showed that the presence of Alfven 
waves in the early Universe leads to specific correlations of the 
CMB in harmonic space: 
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where Q is the power spectrum in the absence of Alfven waves, 
6b and cps are the spherical angles of a PMF direction 5. Here, 



I € f is given by: 



x exp^-2^jv 2 Mkno)Mkno) 



n v +3 



exp | -2— 1 7>(fe7o) ti(folo), 



where 770 is the present conformal time, and denotes the co- 
moving wavenumber of the dissipation scal e, due to photon vis - 
cosity and given by approximately 10/7/i as t (Dur rer et al.|1998| ). 
The damping effec t becomes significant on multipoles i > 500 
( Durrer et al.|1998j ). As shown above, Alfven waves in the early 
Universe produce correlations between harmonic modes sepa- 
rated by A£ = 0, ±2, and Am = 0, ±2. Investigating these im- 
prints, we may impose a constraint on the Alfven waves. In the 
weak Alfven wave limit, the CMB data likelihood can be ex- 
panded as follows: 



d{A v v\) 
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1 d 2 £ 
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(A.5) 



Since all correlations produced by Alfv en waves are propor- 



tional to A v v\, the first term in Eq. A.5 is simply equal to the 



likelihood of the standard cosmological model. The first and sec 
ond derivative of the likelihood are obtained by: 



dA { h dA 2 

where 



= -W 2 ) + CH><<H>, 



(A.6) 



(A.7) 



a is the data vector, consisting of the spherical harmonic coeffi- 
cients, ai m , of the CMB anisotropy data, and C is their covari- 
ance matrix. 

In our analysis, we consider the four foreground cleaned 
CMB maps Commander -Ruler, NILC, SEVEM, and SMICA and 
apply the common mask. We assume the fiducial Planck cosmo- 
logical model, and use MC simulations to estimate the ensemble 
average values for sign al an d noise, as required in Eq. |A.6| The 
quantity C~ l a fro m Eq. A.7 was then deter mined by the messen- 
ger f ield method ([Eisner & Wandelt 2013 ). The CosmoMC pack- 
age ( [Lewis &T3 ridle 2002 ) is then used as a generic sampler 
in order to obtain the posterior probability for the Alfven wave 
parameters {A v v 2 , n v , B , 0s}. 



In Table A.l we show the upper bounds on the Alfven wave 
amplitude A v v z A at various confidence levels, after marginaliz- 
ing over the spectral index n v and the direction 6b, $b- From the 
analysis of the Planck data, we impose an upper bound on the 
Alfven wave amplitude that is tighter than that from the WMAP 
data by more than one order of magnitude. 
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Table A.l. Planck constraints on the Alfven wave amplitude 
A v v\. 



Confidence Level 


68% 




95% 




99.7% 


C-R 


< 0.48 X 10" 


-9 


< 1.01 x 10" 
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< 0.54 x 10" 
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